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These notes introduce and review some of the physical principles underlying the theory of 
astrophysical accretion, emphasizing the central roles of angular momentum transport, 
angular momentum loss, and radiative cooling in determining the structure and evolution 
of accretion flows. Additional topics covered include the effective viscous theory of thin 
disks, classical instabilities of disk structure, the evolution of warped or eccentric disks, 


and the basic properties of waves within disks. 
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I. INTRODUCTION 


Accretion is central to astronomical systems as diverse 
as protoplanetary disks, X-ray binaries, Active Galac- 
tic Nuclei, Gamma-Ray Bursts, and Tidal Disruption 
Events. There are very many important differences be- 
tween these systems, but what they have in common is 
that they all involve accretion, usually from a disk, on to 
a central gravitating object. Because the specific angu- 
lar momentum in a Keplerian potential increases outward 
as h x yr, the physics problem posed by accretion is to 
understand how gas in a rotating flow loses angular mo- 
mentum and moves inward. This problem is generally 
considered to be solved in principle, though important 
aspects remain unclear. The astronomical problem is to 
figure out how the observed manifestations of accretion, 
in the form of resolved images, spectra and time variabil- 
ity, derive from the underlying physics. Many facets of 
this problem are certainly not solved. 


Indispensable references for the student of accretion 
are Jim Pringle’s review of viscous disk theory (Pringle, 
1981), Steve Balbus and John Hawley’s review of turbu- 
lence and angular momentum transport processes (Bal- 
bus and Hawley, 1998), and the textbook by Juhan 
Frank, Andrew King and Derek Raine (Frank et al., 


2002). Gordon Ogilvie’s lecture notes! are extremely 
useful, and I follow his treatment of several classical 
problems here. My goal for these notes is to give an 
accessible introduction that follows a modern point of 
view, in which the key physical processes of angular mo- 
mentum transport and radiative efficiency motivate older 
work on effective viscous disk theory. Also included are 
a smorgasbord of results that are hard to find outside 
the primary literature, extensive (but still incomplete) 
references, and some editorializing on my part. 


ll. ANGULAR MOMENTUM TRANSPORT IN DISK 
ACCRETION 


The angular velocity Q of a particle in a circular orbit 
at distance r from a point mass M in Newtonian gravity 
is, 


M 
where G is the gravitational constant and the subscript 
“K” identifies this as being characteristic of Keplerian 
orbits. The specific angular momentum (i.e. the angular 
momentum per unit mass) is an increasing function of 
orbital radius, 


hx = r°Qg = VGMr. (2) 


This elementary property of orbits in Newtonian grav- 
ity leads immediately to the central problem of accretion 
physics. Very frequently gas that is bound to stars and 
compact objects has specific angular momentum that ex- 
ceeds that of a circular orbit grazing the object’s surface, 
and notwithstanding the complications wrought by gen- 
eral relativity, pressure forces, and so forth, often it orbits 
in a flattened disk-like structure with h ~ hx. An indi- 
vidual parcel of gas must then lose angular momentum 
before it can move to an orbit at smaller r and be ac- 
creted. There are two logical possibilities for how this 
can happen. One possibility is that the gas parcel ez- 
changes angular momentum with another parcel, so that 
one moves in while the other moves out to conserve angu- 
lar momentum. Invariably the viscosity of the gas is too 
small to lead to appreciable angular momentum exchange 
on the time scales inferred for astrophysical systems, so 
to realize this possibility we require that the fluid flow 
exhibits an instability that leads to turbulence and en- 
hanced transport. The other possibility is that the disk 
is not an isolated system and can lose angular momen- 
tum as a whole. This can occur if the disk supports a 
magnetized outflow that exerts a torque back on the disk 
to remove angular momentum. 


l http: //www.damtp.cam.ac.uk/user/gio10/accretion.html 
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FIG. 1 The vertical structure of a thin disk is determined 
by the hydrostatic balance between the vertical component of 
gravity from the central source and a pressure gradient. 


In this Section we will first show that for an impor- 
tant class of disks that are highly flattened, known as 
“thin disks”, the assumption that h ~ hx is easily jus- 
tified. We justify neglecting the regular fluid viscosity 
that results from particle-particle interactions, which is 
negligible in most cases of interest. We then consider the 
linear stability of magnetohydrodynamic (MHD) Keple- 
rian shear flows, along with the case where the disk is 
non-magnetized but self-gravitating. Finally, we will de- 
rive a simple condition for when a thin disk, threaded by 
a large-scale magnetic field, can launch an MHD wind. 


A. Thin disk structure 


Consider a disk of gas that orbits in the z = 0 plane in 
cylindrical polar co-ordinates (r,¢,z). We assume that 
the disk is accreting slowly enough that it is in approxi- 
mate hydrostatic balance, and that the mass of the disk 
is negligible compared to the central mass M. 

To work out the vertical structure of the disk, we look 
at the balance of forces in the vertical direction (shown 
in Figure 1). The vertical component of gravity, 


GM 
3/2 Z, (3) 


z = gsin(ĝ) = —————— 
gz = gsin(0) (42) 


is balanced by a pressure gradient. The vertical struc- 
ture is determined by the equation of vertical hydrostatic 
equilibrium, 


dP 
ae = — pgz. (4) 


In different but reasonable physical circumstances the 
pressure might have important contributions from gas, 
radiation, and magnetic fields. The simplest case is when 
the pressure is dominated by gas pressure, and the verti- 
cal temperature structure is isothermal. This is roughly 
appropriate when the disk is optically thick and heated 
from outside. The equation of state is then, 


P = pc, (5) 


where c,, the sound speed, is not a function of z. If we 
further assume that z < r, the equation of hydrostatic 
equilibrium becomes, 


Ca, -Q2 pz. (6) 


Solving this equation we find that, 


Zz? z? 
(2) = po exp |- k-o- 0 
The density at the disk mid-plane (at z = 0) is pọ, and 
the second equality serves to define the disk scale height, 


Cs 


h= Oe (8) 


It is often convenient to work in terms of the surface 
density 4, defined as, 


L= p(z)dz. (9) 
For the density profile given by equation (7) the relation 
between the mid-plane and surface densities is, 


1s 
Po = me 


The above analysis is justified if h/r = cs/uK < 1, i.e. 
if the disk sound speed is much smaller than the Keple- 
rian orbital velocity. Such disks are described as being 
geometrically thin. 

A geometrically thin disk, by definition, has a mid- 
plane sound speed that is small compared to the Keple- 
rian orbital velocity. In this limit, radial pressure gradi- 
ents do not affect the angular velocity profile of the disk 
at leading order. The radial component of the momen- 
tum equation determines the azimuthal velocity vg of the 
gas, 


(10) 


1dP 
pdr 


2 
US = GM 
2 


(11) 


r ifs 


Approximating dP/dr ~ P/r ~ pc?/r we get, 


1-0(2)]. 02) 


Deviations from Keplerian velocity are thus second or- 
der in h/r, and we can imagine a toy disk model in 
which the angular velocity is Keplerian, the mid-plane 
density and temperature are specified functions of ra- 
dius, and the vertical density profile is gaussian. In the 
right circumstances—when the disk is geometrically thin, 
gas pressure dominated, and isothermal in the vertical 
direction—such a model will be a decent approximation. 

On occasion, one needs a true two-dimensional solu- 
tion for an axisymmetric disk in hydrostatic equilibrium. 


v=o 


Several such solutions are known. Lin et al. (1990) and 
Bate et al. (2002), for example, quote a solution for a disk 
with a polytropic equation of state near the mid-plane, 
and an isothermal atmosphere. This model can be useful 
as a background state when studying wave propagation in 
disks. Fishbone and Moncrief (1976) provide a solution 
for a fluid torus around a black hole. Many numerical 
simulations of black hole accretion are initialized with 
this, or similar, analytic models. 


B. Microphysical viscosity in accretion disks 


To accrete, gas in a Keplerian disk configuration that 
satisfies hydrostatic and rotational equilibrium has to 
lose angular momentum. This can happen in two ways. 
Gas in the disk can lose angular momentum, for example 
when a magnetic field threading the disk exerts a torque 
on the disk surface. Alternatively or additionally, angu- 
lar momentum can be redistributed within the disk, such 
that the inner part of the disk loses angular momentum 
and accretes while the outer part expands to conserve 
angular momentum. 


A rotating fluid redistributes angular momentum due 
to (microscopic) viscosity, but this process is too slow to 
be astrophysically relevant in essentially all disks. For an 
ionized gas of cosmic composition the kinematic viscosity 
at temperature T and density p is (Spitzer, 1962), 


5/2 


T 
v=1.6 x M ar cm? s7}, (13) 
pin 


Here, A is the Coulomb logarithm for proton-proton scat- 
tering. (See Balbus and Henri (2008) for a discussion and 
partial derivation of this result.) It is fiendishly hard 
to measure the density and temperature of most astro- 
physical disks precisely, but there are many observational 
and theoretical ways to get an estimate which is good 
enough for our purposes. For dwarf novae (accreting 
white dwarfs in mass transfer binary systems), for exam- 
ple, we can appeal to state-of-the-art numerical simula- 
tions by Hirose et al. (2014). At a radius of r ~ 10'° cm, 
they obtain characteristic temperatures T = 3 x 104 K 
and densities p ~ 1077 g cm~3. The corresponding vis- 
cosity is v ~ 2500/In A cm? s~!. Anticipating somewhat 
results from §III.C (or on dimensional grounds) we con- 
struct a time scale tT ~ r?/v assuming that the viscosity 
acts diffusively. With these numbers, tT ~ 10° yr. Obser- 
vationally, however, dwarf nova disks are seen to evolve 
dramatically on time scales of just days (for an example 
with nice Kepler data, see Cannizzo et al., 2012). Pretty 
clearly viscosity due to small-scale kinetic processes in 
the gas is not responsible. Similar arguments apply to 
protoplanetary disks and disks in AGN. 


C. The shearing sheet 


The low level of microphysical viscosity in accretion 
disks motivates study of macroscopic instabilities that 
generate turbulence and angular momentum transport. 
The most important of these is the magnetorotational 
instability (MRI; Balbus and Hawley, 1991), which we 
will discuss in §II.D. Our analysis will be based on the 
shearing sheet approximation (Goldreich and Lynden- 
Bell, 1965), which is frequently used in both analytic and 
numerical studies of disks. It is a fluid cousin of Hill’s 
equations (Hill, 1878), developed for the study of lunar 
motion. 

To develop the shearing sheet we start with the inviscid 
momentum equation. In an inertial frame this is just, 


— hy Vy = — - Va: (14) 


The velocity is v, P is the pressure, p is the density, and 
® is the gravitational potential. Equivalently, in terms 
of the Lagrangian or material derivative, 


Soe VG; (15) 


The right hand side expresses the forces acting on a fluid 
element, so transforming to a frame rotating at constant 
angular velocity Q requires adding the Coriolis and cen- 
trifugal forces in the same way as for point mass dynam- 
ics, 
Dv = VP Vo 
Dt p 
Here f is a unit vector in the radial direction. No ap- 
proximation has yet been made. 

Figure 2 shows the geometry of the shearing sheet. The 
idea is to consider a local patch of the disk, centered at ro, 
that co-rotates with the fluid at angular velocity Qo. The 
dynamics is modeled in a Cartesian co-ordinate system, 
within which disk quantities are replaced by constants or 
via first-order Taylor expansion about ro. Mathemati- 
cally, assume that the angular velocity can be approxi- 
mated as a power-law, 


Q(r) = Qo (E) (17) 


To 


20 x v +r t. (16) 


q = 3/2 corresponds to the Keplerian case (equation 1). 
We define a co-rotating co-ordinate system, 


z=r-— To, (18) 
y = rb —Qot, (19) 
and sum up the radial component of the gravitational 


and centrifugal forces, keeping only the first order term 
in z/r9, 


—q 
-VP + ri = -rog (=) +r, (20) 


re WqQN2er. (21) 


Adding in the vertical gravitational acceleration, with a 
value appropriate to the center of the shearing sheet (us- 
ing equation 6), the momentum equation in the shearing 
sheet approximation is, 


D P 
om = = — 2 x v + 2gMxk — 022%. (22) 


Å and Z@ are unit vectors in the x (radial) and z (vertical) 
directions respectively. This three-dimensional version is 
called a shearing box. 

In addition to analytic applications, the shearing box is 
a useful tool for local numerical simulations of accretion 
disks. The key point is that because the angular veloc- 
ity is constant across the box, the dynamical time scale 
Qt is not a function of the radial co-ordinate x. This 
makes it possible to construct self-consistent shearing pe- 
riodic radial boundary conditions. The construction is 
illustrated in Figure 3 (Hawley et al., 1995). The active 
computational domain is bordered radially by copies of 
itself, which shift azimuthally over time at a rate that 
reflects the shear across the domain. For a domain of 
radial extent Lz, the radial boundary conditions can be 
written as (Hawley et al., 1995), 


f(a, y, 2) = f(x + Le, y — WoLct, z), (23) 
Vy (2, Y, 2) = Vy (x as Lr, y= qQoLzt, z) F qQoLz, (24) 


where the first expression applies to all flow variables 
apart from the azimuthal velocity vy, which requires a 
boost to account for the shear across the box. 

The shearing box is a powerful tool. For analytic cal- 
culations, it captures much of the essential small-scale 
dynamics of disks, while being easier to deal with than 
global equations. For numerical work, on the other hand, 
equation 22 is not much easier to solve than its global 
counterpart. The advantage for computation studies de- 
rives mainly from the ability to use shearing periodic 
radial boundary conditions, which are (relatively) non- 
coercive. If, instead, one attempts to simulate a small 
cylindrical patch of disk, it proves hard to devise bound- 
ary conditions which do not change the flow dynamics. 
The shearing box can be readily adapted to include mag- 
netic fields (we will do so shortly), an energy equation, 
and radiation fields (Hirose et al., 2009; Jiang et al., 
2013). It can also be extended to model disks whose back- 
ground states are eccentric (Ogilvie and Barker, 2014) or 
warped (Paris and Ogilvie, 2018). 

The shearing box approximation can also set traps for 
the unwary. At the most basic level, although we defined 
x = r — ro, going to the shearing sheet introduces a sym- 
metry between +a and —z, such that there is no physical 
sense in which one direction is “inward” and the other 
“outward”. This means, for example, that while an MRI 
simulation may develop sustained turbulence and mea- 
surable stresses, the gas does not actually move radially. 
More perniciously, this symmetry means that a shearing 
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FIG. 2 Geometry of the local or shearing sheet model for accretion disks. The shearing sheet is a Cartesian (x, y) representation 
of a small patch of the disk, centered at r = ro, that has an angular velocity Q = Qo. The disk shear is approximated by a 


linear function across the sheet. 


box threaded by a net vertical magnetic field can develop 
a wind solution in which the gas flows to +a (say) above 
the mid-plane, and to —x below the mid-plane. This sort 
of solution is unphysical in a global context, but there is 
nothing to prevent it developing within a shearing box 
(e.g. Bai and Stone, 2013). Another difficulty comes from 
attempts to model radial gradients in disk properties, or a 
vertical gradient in Q. These gradients are physically im- 
portant for various disk instabilities, but they are hard to 
incorporate consistently within the shearing box formal- 
ism (Latter and Papaloizou, 2017; McNally and Pessah, 
2015). 


D. Magnetorotational instability 


The magnetorotational instability (MRI) is the most 
important process that can generate turbulence and an- 
gular momentum transport in accretion disks (Balbus 
and Hawley, 1991, 1998). The MRI is a local, linear in- 
stability, that is present in accretion disks provided that, 


(i) There is differential rotation, with dQ/dr < 0. Ex- 
cept in boundary layers, all disks satisfy this con- 
dition. 


(ii) “Weak” magnetic fields are present. The MRI can 
be damped by microphysical viscosity for very weak 
fields (too weak to be relevant to any disks apart 
perhaps from those around primordial stars; Kro- 
lik and Zweibel, 2006), and stabilized (though not 
necessarily entitrely shutdown) by magnetic ten- 
sion effects in disks where the magnetic pressure is 
larger than the thermal pressure (Das et al., 2018; 
Pessah and Psaltis, 2005). 


(iii) The ionization fraction is high enough to couple the 
magnetic field to the fluid. This is the most impor- 
tant caveat. Although very low ionization fractions 
suffice, protoplanetary disks can be so neutral as to 
call into question the importance of the MRI (Gam- 
mie, 1996). 


These prerequisites are weak, and most astrophysical 
disks satisfy them. Accordingly, we will skip over the 
tricky and rich question of whether non-magnetic non- 
linear? fluid instabilities would exist in notional non- 
magnetic disks, and proceed directly to the MHD case. 
Lyra and Umurhan (2019) review potential fluid instabil- 
ities in protoplanetary disk systems where the viability 
of the MRI as a transport mechanism remains doubtful. 

The physical origin of the MRI is shown in Figure 4, 
for the conceptually simplest case of a disk threaded by a 
weak, initially vertical magnetic field. Consider a radial 
perturbation to the magnetic field—exaggerated in the 
cartoon—that results in linkage between an inner (blue) 
fluid element and an outer (green) one. As the disk ro- 
tates, differential rotation causes the fluid elements, ini- 
tially at the same azimuth, to become azimuthally sep- 
arated. This separation is opposed by magnetic tension, 
which leads to a force that acts to decrease the angu- 
lar momentum of the inner fluid element and increase 
that of the outer element. Because angular momentum 
is an increasing function of radius in Keplerian disks, this 


2 Keplerian flows are linearly stable by the Rayleigh criterion, as 
the specific angular momentum is an increasing function of ra- 
dius. 
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FIG. 3 Illustration of how radial boundary conditions are imposed in computational studies that use the shearing box approx- 
imation (Hawley et al., 1995). Conceptually, we imagine that the simulation domain is bordered radially by identical copies of 
itself, which shear azimuthally over time. Practically, these boundary conditions are often implemented using time-dependent 
copying of the fluid quantities in the active zones into one or more layers of “ghost” zones, interpolating as necessary when the 


required azimuthal shift is not an integer number of grid cells. 


transfer of angular momentum causes the fluid elements 
to separate further radially, signalling an instability. 

The existence of an instability in magnetized rotating 
flows was demonstrated by Velikhov (1959) and by Chan- 
drasekhar (1961), but the importance of this instability 
for accretion disks was not recognized. Safronov (1972) 
came close, but no cigar. The MRI was rediscovered and 
applied to accretion disks in breakthrough work by Bal- 
bus and Hawley (1991). 


1. The MRI dispersion relation 


The dispersion relation for the MRI can be derived in 
several distinct ways. Here, following Fromang (2013), 
we work through a rather explicit calculation of a differ- 
entially rotating disk containing a purely vertical mag- 
netic field, whose stability we assess in the shearing sheet 
approximation. The disk has a power-law angular veloc- 
ity profile, Q œ r74, and a uniform vertical magnetic 
field Bo. Radial and vertical variations in density, and 
the vertical component of gravity, are ignored. The equa- 
tion of state is isothermal, P = pc?, with c, a constant. 
As with any linear stability analysis the goal is to find 
out whether infinitesimal perturbations to the equilib- 


rium are stable, or whether instead they exhibit expo- 
nentially growth, implying instability. 

The relevant equations are the continuity equation, the 
induction equation in the ideal magnetohydrodynamic 
(MHD) limit, and the momentum equation (including 
MHD forces) in the shearing box approximation, 


Op 
az tV (ev) =9, (25) 
Ov 1 1 
290 X v + 2qN2rx, (26) 
oB 
op TV x xB). (27) 


The equilibrium state has uniform density, p = po, and 
uniform magnetic field B = (0,0, Bo). In this setup there 
are no pressure or magnetic forces, and the initial velocity 
field is set by the balance of the Coriolis and centrifugal 
terms, 


2 X v = QNR. (28) 
The initial velocity field is then, 


v = (0, —qQox, 0). (29) 
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FIG. 4 The operation of the magnetorotational instability (MRI) for a perturbed vertical magnetic field in ideal MHD. The 
MRI is a fluid instability, but the dynamics is closely related (and in some cases mathematically equivalent) to that of two 
orbiting point masses connected by a spring (Balbus and Hawley, 1992). 


As expected given that we are working in the shearing 
sheet, there is a linear shear around radius rọ. 

To determine the stability of this system, we write the 
fluid variables as the sum of their equilibrium values plus 
perturbations. For the MRI, it suffices to consider a per- 
turbation which depends only on z and t. For the veloc- 
ity, for example, we take, 


Uz = vi,(z, t), (30) 
Vy = —qox F vy(z,t), (31) 
vz = v,(z, 2), (32) 


and do similarly for the magnetic field and density. These 
expressions are substituted into the continuity, momen- 
tum and induction equations, discarding any terms that 
are quadratic in the primed quantities. This procedure 
leads to seven equations, but to derive the MRI we need 
only a subset of them, namely the x and y components 
of the momentum and induction equations. The four lin- 
earized equations are, 


Ov! Bo OB' 
x = ng 2Q 7 
Ot Ampo Oz 09y» 3] 
ðv! B) OB 
y L aov, = 2 o! 4 
Ot q OUy Ar po Oz 0Vz> (3 ) 
OB! ðv! 
zZ = Bao 35 
at 0 Az , ( ) 
OB, Ov, ; 
Begg ea (35) 


These linearized differential equations are then converted 
into algebraic equations by taking the perturbations to 


have the form, 
B, = B elik), (37) 


where w is the frequency of a perturbation whose vertical 
wave-number is k. The time derivatives give us factors 
of iw, while the spatial derivatives give us ik. The four 
algebraic equations we end up with are, 


, BoB, 

iwv, = —ik ten + 2Dov,, (38) 
_ BoB, 

iwu, = —tk Ta + (q—2)Qov},, (39) 
iwB! = —ikBov',, (40) 


(Here we’ve dropped bars on the variables.) The final 
step is to eliminate the perturbation variables between 
the equations, leaving us with the dispersion relation for 
the MRI. It takes the form, 


wt [2k?v% + 2(2 — q)92] + 
k?v [k v3 — 2492] = 0, (42) 


where v4, = Bê /(4rpo) is the Alfvén speed in the initial 
state. 

The MRI dispersion relation is a quadratic in w?. It 
is plotted as a function of kva, for different values of 
the shear parameter q, in Figure 5. For the q = 3/2 
case that characterizes Keplerian disks w? < 0, implying 
instability, for all scales larger than some minimum (i.e. 
for kva smaller than some maximum). Graphically, it 
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FIG. 5 The MRI dispersion relation in ideal MHD is plotted 
for different rotation laws. The lower red curve shows the 
q = 3/2 case appropriate for a Keplerian disk. The blue 
curves are plotted, from top downwards, for q = —1/2, q = 0, 
q = 1/2, and q = 1. Instability (w?/§ < 0) is present at 
sufficiently small k for q > 0. 


can be seen that the instability disappears for q < 0. 
Requiring that w? < 0 the instability criterion is, 


Q2 
<0. (43) 
r 


In the limit where the vertical field By —> 0, the condition 
for instability is just that the angular velocity decreases 
outward. This is distinct from the Rayleigh stability cri- 
terion for a strictly unmagnetized fluid, which instead 
requires that the specific angular momentum r?Q. of the 
fluid decrease outward for an instability (see e.g. Pringle 
and King, 2007). 

A number of other important quantities can be 
straightforwardly derived from equation (42). Solving 
for dw? /d(kva) = 0 on the unstable branch of the disper- 
sion relation yields the scale with the fastest MRI linear 
growth rate. For a Keplerian disk the answer is, 


V15 


(ku A) max = a o- (44) 
The growth rate at this spatial scale is, 
3 
önal = g2% (45) 


This is a fast growth rate indeed! The most unstable 
MRI modes grow by a factor exp(37/2) ~ 10? per orbit. 

Returning to equation (42) we can set w? = 0 and 
find kerit, the largest wavenumber (smallest spatial scale) 
which is unstable. Specializing again to the Keplerian 
case we find, 


keritva = V390. (46) 


This result can be used to approximately quantify our 
earlier assertion that the MRI can be stabilized if the 


vertical field is too strong. For a strong field, the small- 
est unstable scale Ari, may exceed the vertical thickness 
of the disk, leading to a situation where the MRI is sta- 
bilized on account of the disk being too thin to support 
any unstable modes at all. For an estimate, we can set, 


20 


crit a? 
kerit 


= 2h, (47) 
which leads to an estimate of the stabilizing field as, 


2 
DG max = — PCs. (48) 
T 
Defining the plasma 8 parameter to be the ratio of the 
gas pressure to the magnetic pressure, 


8r P 


= 49 
B= or (49) 
and making use of h = c,/Qo, we can re-express equa- 


tion (48) as, 


Qn? 

pem (50) 
Having ignored vertical density gradients—which by def- 
inition are significant on scales of the disk scale height 
h—this can only be an estimate. With that caveat, 
the conclusion is that disks are unstable to the MRI 
provided that the vertical magnetic field is significantly 
sub-thermal, in the sense of the magnetic pressure being 
smaller than the thermal pressure. 

Somewhat surprisingly, relaxing the severe simplifica- 
tions that we have made in the above analysis—a local 
shearing sheet, with a simple field geometry, simple per- 
turbations, and no consideration of the flow energetics— 
does not materially change the answer. The review by 
Balbus and Hawley (1998), and more recent work by Lat- 
ter et al. (2015), are good places to start for more com- 
plete analyses of the MRI. 

The linear analysis says nothing about whether the 
MRI leads to physically significant levels of turbulence 
and angular momentum transport in accretion disks. For 
this, we must turn to simulations, in either shearing box 
(Brandenburg et al., 1995; Hawley et al., 1995) or global 
(Armitage, 1998) geometries. A local simulation of the 
MRI for a vertically stratified isothermal disk (Simon 
et al., 2012) is shown in Figure 6. The most basic sim- 
ulation results, known now for more than twenty years, 
include the fact that the MRI can act as a dynamo, sus- 
taining turbulence and magnetic fields against dissipa- 
tion in a domain with no external currents. Angular 
momentum is transported primarily due to the action of 
MHD (Maxwell) stresses. Stronger levels of turbulence 
and transport occur if the disk is threaded by a net ver- 
tical magnetic field. There remain open questions, of un- 
certain physical significance, regarding the convergence 
of the MRI in local domains with simple physics (at the 
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FIG. 6 Rendering of the current density from a local shearing 
box simulation of MRI disk turbulence. Based on simulations 
from Simon et al. (2012). 


least, very high resolution is needed; Ryan et al., 2017), 
but the bulk of recent work has turned to models with 
more complete treatments of thermal, radiative, relativis- 
tic, and plasma physics. Jiang et al. (2019a), for exam- 
ple, use radiation magnetohydrodynamic simulations to 
study the global evolution of disks around supermassive 
black holes. 


E. Self-gravity 


In addition to the MRI, there are other linear insta- 
bilities of accretion disk flow. Most of these either occur 
under less general circumstances than the MRI, or have 
much lower growth rates, making them less important. 
Self-gravity is a partial exception. The dynamics of a 
fluid disk that is massive and / or cold is strongly modi- 
fied by self-gravity if the Toomre Q parameter, 


<1. (51) 


Here « is the epicyclic frequency, K? = 40? + 2rQdQ/dr. 
A similar result (and the eponymous work; Toomre, 1964) 
applies to particle disks, for which the sound speed in the 
equation should be replaced by the velocity dispersion. 
Conditions where Q ~ 1 are likely to occur at sufficiently 
large radii in Active Galactic Nuclei (AGN) disks (Shlos- 
man et al., 1990), and during star formation. Gravita- 
tional instability can in turn lead to angular momentum 
transport, or fragmentation (for a review, see Kratter 
and Lodato, 2016). 

We can gain intuition into where Q comes from using 
a simple time scale argument. Pressure will prevent the 
collapse of a patch of the disk, with scale Ar, when the 
sound-crossing time Ar/c, is small compared to the free- 


fall time ,/Ar?/GAr?d. (Ignoring factors of the order 


of unity.) Equating these time scales the minimum scale 
of collapse can be estimated as Ar ~ c?/GX. On larger 
scales, collapse can be stopped by shear if the free-fall 
time is long compared to the time scale for radial shear 
to separate neighboring fluid elements. In a Keplerian 
disk the shear time scale is ~ Q7!. For a disk that is 
marginally unstable the minimum scale set by pressure 
equals the maximum scale set by shear. This condition 
implies that marginal stability occurs when c,Q/GE ~ 1, 
as quoted above. 


1. Dispersion relation 


Proceeding more formally (Pringle and King, 2007) we 
consider a razor-thin circular gas disk with uniform sur- 
face density Xo and sound speed cs in the z = 0 plane. 
In cylindrical polar coordinates (r, œ, z) the density of the 
disk is given by, 


Po (r; d, z) = X06 (z), (52) 


where 6(z) is a Dirac delta-function. The velocity field 
is, 


vo(r, ¢, z) = (0, rQ, 0). (53) 


The angular velocity Q(r) is not required to be Keplerian, 
but for circular orbits the centrifugal force must balance 
gravity. If the gravitational potential is Bo, 


Pr= (54) 


The assumptions of initially constant density and sound 
speed mean that pressure gradient forces do not enter 
into the calculation of the equilibrium state. 

The stability analysis proceeds in a similar way as for 
the MRI, except that we now ignore magnetic fields while 
including both the gravity of the central object and the 
self-gravity of the disk itself. The equations are the conti- 
nuity and momentum equations, together with Poisson’s 
equation for the gravitational field, 


oD») 

OE +V- (Xv) =0, (55) 

Ov _ Vp 

at Vy =-F- VS, (56) 
V2 = 47Gp. (57) 


(58) 


defined in terms of the pressure p and surface density © 
in the usual way. 


We consider infinitesimal axisymmetric perturbations 
to the equilibrium state, 


£= Ep + Vale, t), (59) 
P= Po + pir, t), (60) 
d= Do Ie (r, t), (61) 
v = vo + [vr (r, t), dvg(r, t), OJ, (62) 


that have a spatial and temporal dependence given by 
(using the surface density as an example), 


£ (r,t) x expli(kr — wt)]. (63) 


Here k is the spatial wavenumber of the perturbation 
(related to the wavelength via A = 27/k) and w is the 
temporal frequency. Making one further approximation, 
we assume that for the perturbations of interest, 


kr > 1. (64) 


This amounts to considering disturbances that are small 
compared to the radial extent of the disk. 

We now substitute the expressions for the surface den- 
sity, pressure, gravitational potential and velocity into 
the fluid equations, discarding any terms we encounter 
that are quadratic in the perturbed quantities. For the 
continuity equation this yields, 


1 
sey +v, Xo ( + it) =0, (65) 
Tr 


which simplifies further in the local limit (kr >> 1) to 
—wd4 + kv, dig = 0. (66) 


Deriving the analogous algebraic equations from the mo- 
mentum equation requires us to express the convective 
operator (v - V)v in cylindrical coordinates. This takes 
the form, 


vr _ U¢ Ou, ðv, v5 
(v Viv = Ẹ a | r ð$ T Uz əz ee 
Ove Zi Vo Ove Ove UrVe 
or Op r Og *Oz r” 
Ov, ve Ov; : Ov, 
Ur Or A ad F UZ d (67) 


ii, — svg = ae (68) 
0 


iwdvg + v, Ẹ ! =o] =0, (69) 


where the two equations come from the radial and az- 
imuthal components respectively. 

The next step is to relate the perturbations in pres- 
sure and gravitational potential expressed on the right- 
hand-side of equation (68) to perturbations in the surface 
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density. For the pressure term this is straightforward. 
Equation (58) implies that, 


cikd). (70) 
0 


Dealing with the potential perturbations requires more 
work. Starting from the linearized Poisson equation, 


V? = 4rGd16(z), (71) 


we write out the Laplacian explicitly and simplify making 
use of the fact that for short wavelength perturbations 
kr >> 1. This yields a relation between the density and 
potential fluctuations, 

d?’ 

dz2 
For z Æ 0 the only solution to this equation that remains 
finite for large |z| has the form, 


®, = Cexp[—|kzl], (73) 


where C remains to be determined. To do so we integrate 
the Poisson equation vertically between z = —e and z = 
+e, 


k?@, + 4nGd16(z). (72) 


+e +e 
J V?O,dz = / 4nG16(z)dz. (74) 
Noting that both 0?®, /Ox? and 07, /Oy? are continuous 
at z = 0, whereas 07/02? is not, we obtain, 


aby 
dz |e 


Taking the limit € > 0 we find that C = —27G¥1/|k|, 
and hence that the general relation between potential and 
surface density fluctuations on the z = 0 plane is, 


27Gd, 
@,=—- : 


+e 
= 4rG¥;. (75) 


Taking the radial derivative, 


d®, 2rikGX: 
which allows us to eliminate the potential from the right- 
hand-side of equation (68) in favor of the surface density. 
The result is, 


Dz: _ 2rikGX 

SA iky 4 k (78) 
Finally, we are ready to derive the functional relationship 
between w and k (the dispersion relation). Eliminating 
v, and dug between equations (66), (69), and (78) we find 
that, 


iw, — 20dv¢g = 


w? = k? + Êk? — 2nGDolkl, (79) 
where the epicyclic frequency « is defined as, 
dQ 
2 = 40? + 270. 
K + 2r Tr (80) 


In a Keplerian potential k? = Q2. 
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FIG. 7 The dispersion relation for a self-gravitating accretion 
disk in the Keplerian limit is plotted for different values of 
the Toomre Q parameter, Q = 0.75, Q = 1 and Q = 1.5. 
Instability is present at some scales for Q < 1. The blue dot- 
dashed lines show the contributions from pressure (dominant 
at large |k|), rotation (dominant as |k| — 0, and self-gravity 
(destabilizing). 


2. Toomre Q 


Let us consider, for simplicity, a Keplerian disk. Not- 
ing that h = c,/Q, we can write the dispersion relation 
in the form, 


w? 2h|k| 325 

where Q = c,Q/(7G™) is the dimensionless Toomre Q 
parameter that we previously deduced using a time scale 
argument. It is plotted in Figure 7 for different values 
of Q. Each of the terms on the right-hand-side has a 
simple physical interpretation. The constant term de- 
scribes the effect of rotation, which stabilizes all scales 
but which is dominant at small hk] (i.e. at large spatial 
scales). The quadratic term describes the effect of pres- 
sure, which has the opposite tendency and preferentially 
stabilizes small spatial scales. The linear term, which 
describes self-gravity, is negative and thus destabilizing. 
The strength of self-gravity is fully parameterized by the 
value Q, with small enough Q leading to w? < 0 and un- 
stable modes. It is easy to verify that the criterion for 
instability corresponds to Q < 1. 


3. Angular momentum transport or fragmentation 


The onset of self-gravity in an accretion disk can lead 
to angular momentum transport or to fragmentation. 
Loosely speaking, fragmentation is the outcome for disks 
that either cool too quickly (Gammie, 2001; Rice et al., 
2005), or that are fed mass from infall at too high a rate 
(Kratter et al., 2010). Kratter and Lodato (2016) dis- 
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cuss the quantitative thresholds, determined from simu- 
lations, for these outcomes. 


F. MHD disk winds 


Several physical processes can lead to mass loss in a 
disk wind from geometrically thin accretion disks, includ- 


ing, 


e A vertical gradient of thermal pressure. A suffi- 
ciently strong thermal pressure can develop, even in 
a thin disk, if high energy radiation heats the disk’s 
surfaces to a temperature where cs ~ vg. The hot 
gas then has positive total energy and can escape 
to infinity. Thermal or photoevaporative winds may 
be important in AGN (Begelman et al., 1983) and 
in protoplanetary disks (Bally and Scoville, 1982), 
where they contribute to the dispersal of the disk at 
the end of the protoplanetary disk phase (Alexan- 
der et al., 2014; Ercolano and Pascucci, 2017). 


e Radiation pressure acting on spectral lines. The 
physics of line-driven disk winds is an extension 
of the accepted theory for mass loss from mas- 
sive stars (Castor et al., 1975). For accretion 
disks, line-driving is efficient for AGN disks that 
are strong emitters of ultraviolet radiation (Proga 
et al., 2000). (Radiation absorbed by the con- 
tinuum can also be important, but typically only 
under conditions where the disk is geometrically 
thick.) 


e MHD acceleration, which in general involves con- 
tributions from centrifugal acceleration (Blandford 
and Payne, 1982) along poloidal magnetic field 
lines (i.e. in the (r,@) plane of spherical polar co- 
ordinates), and from a gradient of toroidal mag- 
netic pressure (Lynden-Bell, 1996). 


In addition to these processes, which in principle can 
drive a wind from a broad range of disk radii, there are 
others whose applicability is limited to the innermost disk 
region. Disk magnetic fields interacting with a spinning 
black hole can extract spin energy via the Blandford- 
Znajek effect (Blandford and Znajek, 1977), while for 
accreting objects with a material surface the interaction 
between the inner disk and a magnetosphere can generate 
outflow (Shu et al., 1994). Inner disk processes are impli- 
cated in the formation of well-collimated jets in accreting 
systems. 

Of these mass loss processes, magneto-centrifugal disk 
winds are particularly interesting because they can ex- 
tract angular momentum as well as mass from the un- 
derlying disk flow. The basic idea is that gas, accelerat- 
ing outwards along a poloidal field line, rotates with the 
angular velocity of the gas at the field line’s foot point 


in the disk. The departing gas thus gains angular mo- 
mentum, which is removed from the disk by a magnetic 
torque that can be considered to act on the disk’s sur- 
face. The angular momentum loss leads to inflow of gas 
in the disk, independent of internal angular momentum 
transport processes such as the MRI or self-gravity. 

The clearest description of the physics of MHD disk 
winds that I’m aware of can be found in Spruit (1996). 
Here, we content ourselves with a derivation of one of the 
basic properties of a Blandford and Payne (1982) disk 
wind—the existence of a critical inclination angle for a 
poloidal disk-threading field to launch a cold wind. 


1. Blandford-Payne disk winds 


The geometry of a Blandford-Payne wind (Blandford 
and Payne, 1982) wind is illustrated in Figure 8. We en- 
visage a Keplerian disk threaded by a large-scale poloidal 
magnetic field, in the limit of ideal MHD. Within the disk 
the energy density in the magnetic field, B?/87, is usu- 
ally smaller than pc?, the thermal energy. Due to flux 
conservation, however, the energy in the vertical field 
component, B?/87, is roughly constant with height for 
z < r, while the gas pressure decreases rapidly (for a 
thin isothermal disk, as a gaussian with a scale height 
h <r). This leads to a region above the disk surface 
where magnetic forces dominate. The magnetic force per 
unit volume can be written as the sum of a magnetic 
pressure gradient and a force due to magnetic tension, 


JxB B? B-VB 
= 2 
c a (=) Eo (ay 
where the current, 
c 
J= —VxB. (83) 
AT 


In the disk wind region where magnetic forces dominate, 
the requirement that they exert a finite acceleration on 
the low density gas can only be satisfied if the force ap- 
proximately vanishes, i.e. that, 


JxBw0. (84) 


The structure of the magnetic field in the magnetically 
dominated region is then described as being “force-free” , 
and in the disk wind case (where B changes slowly with 
z) the field lines must be approximately straight to en- 
sure that the magnetic tension term is small. If the field 
lines support a wind, the force-free structure persists up 
to where the kinetic energy density in the wind, pv”, first 
exceeds the magnetic energy density. This criterion de- 
fines the Alfvén surface. Beyond the Alfvén surface, the 
inertia of the gas in the wind is sufficient to bend the 
field lines, which wrap up into a spiral structure as the 
disk below them rotates. 
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Magneto-centrifugal driving can launch a wind from 
the surface of a cold gas disk if the magnetic field lines 
are sufficiently inclined to the disk normal. The criti- 
cal inclination angle in ideal MHD can be derived via an 
exact mechanical analogy. To proceed, we note that in 
the force-free region the magnetic field lines are (i) basi- 
cally straight lines, and (ii) enforce rigid rotation out to 
the Alfvén surface at an angular velocity equal to that 
of the disk at the field line’s footpoint. The geometry is 
shown in Figure 8. We consider a field line that inter- 
sects the disk at radius r9, where the angular velocity is 
Qo = \/GM/r@, and that makes an angle 8 to the disk 
normal. We define the spherical polar radius r, the cylin- 
drical polar radius w, and measure the distance along the 
field line from its intersection with the disk at z = 0 as s. 
In the frame co-rotating with Qo there are no magnetic 
forces along the field line to affect the acceleration of a 
wind; the sole role of the magnetic field is to constrain 
the gas to move along a straight line at constant angular 
velocity. Following this line of argument, the condition 
for the acceleration of the wind can be described in terms 
of an effective potential, 


GM 
r(s) 
that is the sum of the gravitational potential and the 
centrifugal potential in the rotating frame. 
Written out explicitly the effective potential is, 
GM 
(s2 + 2sro sin 0 + r2)1/2 


Bels) = -2A — 5030s), (85) 


Der (s) = 
= 1 Q2 (r + inĝ 2 86 
5) 0 (ro ssin ) * ( ) 


This function is plotted in Figure 9 for various values 
of the angle 0. If we consider first a vertical field line 
(9 = 0) the effective potential is a monotonically in- 
creasing function of distance s. For modest values of 
0 there is a potential barrier defined by a maximum at 
some S = Smax; While for large enough 0 the potential 
decreases monotonically from s = 0. In this last case 
purely magneto-centrifugal forces suffice to accelerate a 
wind off the disk surface, even in the absence of any ther- 
mal effects. We compute the critical inclination angle 
Orit, defined as the minimum angle that allows magneto- 
centrifugal wind driving, via the condition, 


3? Dog 
Os? 


=0. (87) 
s=0 


Evaluating this condition, we find, 
1 — 4sin? Orit = 0 
=> berit = 30°, (88) 


as the minimum inclination angle from the vertical 
needed for unimpeded wind launching in ideal MHD. 
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FIG. 8 Illustration (adapted from the pedagogical review by Spruit, 1996) of a magnetized disk wind driven by centrifugal 
force. Left panel: the region above the disk surface is force-free. If magnetic field lines threading the disk are inclined by a 
large enough angle with respect to the vertical, centrifugal force can accelerate gas along the field lines even in the absence of 
a pressure gradient. The dynamics is equivalent to a mechanical system of a bead (the gas) on a rigid wire (the magnetic field 
line) that rotates with the angular velocity of the disk at the foot point. Right panel: the geometry for calculating the critical 


inclination angle for cold magneto-centrifugal wind launching. 
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FIG. 9 The effective potential for a cold magnetocentrifugally 
driven disk wind. Acceleration occurs for field lines that are 
inclined by 30° or more away from the vertical. 


The rigid rotation of the field lines interior to the 
Alfvén surface means that gas being accelerated along 
them increases its specific angular momentum. The mag- 
netic field, in turn, applies a torque to the disk that re- 
moves a corresponding amount of angular momentum. If 
a field line, anchored to the disk at radius 79, crosses the 
Alfvén surface at (cylindrical) radius r4, it follows that 
the angular momentum flux is, 

y = Mor, (89) 
where M,, is the mass loss rate in the wind. Removing 
angular momentum at this rate from the disk results in 
a local accretion rate M = 2L,,/Qorz. The ratio of the 


disk accretion rate to the wind loss rate is, 
(90) 


If ra substantially exceeds ro (by a factor of a few, which 
is reasonable for detailed disk wind solutions) a relatively 
weak wind can carry away enough angular momentum to 
support a much larger accretion rate. 


Il. EFFECTIVE VISCOUS THEORY OF ACCRETION 
DISKS 


The MRI, self-gravity, and disk winds are physical pro- 
cesses that lead to angular momentum transport and evo- 
lution of accretion disks. In many cases, it is not possible 
to simulate these processes with sufficient fidelity, or over 
enough radial range and for long enough, to compare 
against observations. It can therefore be necessary to 
turn to viscous disk theory, historically developed much 
earlier, which can be used to model long-term disk evolu- 
tion (Lynden-Bell and Pringle, 1974; Shakura and Sun- 
yaev, 1973). 

The basic assumption of an effective viscous disk the- 
ory is that angular momentum transport within the disk 
can be represented approximately as a fluid viscosity, 
using some parameterization that does not explicitly in- 
volve B. The continuity and momentum equations are 
then, 

ĉr +V-(pv) =0, 
OY pgp AVG Ve 
ðt p p 


(91) 


(92) 


Here p is density, v is velocity, P is pressure, ® is gravita- 
tional potential, and T is stress (represented by a tensor). 
We further assume that we are dealing with a geomet- 
rically thin disk for which the pressure gradient term is 
small. We can then make progress using simplified ver- 
sions of the equations for mass and momentum conser- 
vation. 


A. One dimensional time-dependent disk evolution 


The fluid equations (92) apply generally. We first spe- 
cialize to the case of a geometrically thin, circular, and 
planar disk, and derive a one-dimensional (in radius r) 
evolution equation. We then look for time-independent 
and time-dependent solutions. 


1. 1D evolution equation 
In cylindrical polar co-ordinates (r, ¢, z) the continuity 
equation is, 


dp 19 1a ð 
BE ra e ag OAT a ee O88) 


Integrating over ¢ [0,27] and over z [—oo 
term becomes, 


ð 2T oo ð 
I f F pdzdod = a (275) , (94) 


where È(r, t) is the azimuthally averaged surface density. 
We are working toward an evolution equation for this 
quantity. Integrating the second term, 


2T 
-f rpvrdzdġ = i (95) 
r = 


where F, the radial mass flux, 


oo] the first 


27 co 
FH I rpu,dzdd = 2n Dū, (96) 
0 —co 


can be written in terms of the surface density and the 
density-weighted radial velocity v,. On integration, the 
third and fourth terms vanish (the latter assuming that 
there is no mass loss from the surfaces of the disk) leav- 
ing, 
+ ie ema) =0, (97) 
as the one-dimensional version of the continuity equation. 
Dealing with the momentum equation in the same way, 
we can simplify the algebra by assuming at the outset 
that the disk is axisymmetric such that the orbital veloc- 
ity, 


ve =7Q, (98) 
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depends only on the angular velocity Q = Q(r) of circu- 
lar orbits in a fixed gravitational potential ®. The only 
non-zero terms in the ¢ component of the momentum 
equation then come from v: Vv and V-T. Looking up 
the forms for these in cylindrical polar co-ordinates the 
surviving terms are, 


dvg | Ugur 
v- Vv]; = vr Pa (99) 
1 1 3 1 ôTyz 
S E Cea re 1 
y o põr (Trg) + Oz (100) 


Multiplying these expressions through by rp the az- 
imuthal component of the momentum equation takes the 
form, 


(101) 


where the specific angular momentum h is defined as, 


h= rQ. (102) 
Finally we multiply equation (101) through again by r, 
and integrate over ġ and z. If Ty, vanishes as z —> +00 
the result is, 


dh OG 
Sea ees 1 
dr Or’ (103) 
where F is given by equation (96) and, 
27 o0 
=- | l r’T,.gdzdd, (104) 
0 —oo 


is the viscous torque. 

Getting to this point from the conservation laws ex- 
pressed in equation (92) requires only some rather trans- 
parent assumptions: axisymmetry, a time-independent 
potential, and no mass or angular momentum loss from 
the disk surfaces. One could stop there, and consider T;.¢ 
to be the key quantity whose dependence on disk condi- 
tions needs to be determined. Conventionally, however, 
we instead write the torque in terms of an effective fluid 
viscosity that follows a Navier-Stokes form. For a fluid 
with viscosity u and bulk viscosity up we have, 


2 
T=p [vv + (Vv)"| + (1 — 3H) (V-v)I. (105) 
To leading order the divergence of a thin disk velocity 
field vanishes, so we don’t have to worry about bulk vis- 
cosity at all. The rø component of the stress is, 


dQ 


pe = ur, 1 
Trp = ur (106) 


Defining the kinematic viscosity (later just “the viscos- 
ity”) v as, 


2T 
ay ie pdzdd, (107) 


the viscous torque has a fairly intuitive form that is the 
product of the circumference, the viscous force per unit 
length, and the lever arm, 


Q 
G=-—2zr- a r. 


Tr (108) 
Equation (103) is then, 
dh ld 3 dQ 
ey dr rdr (>r dr ) (uos) 


Given the aforementioned assumptions, this equation ex- 
presses angular momentum conservation for a viscous 
fluid in a disk geometry. 

Eliminating the radial velocity 0, between equa- 
tion (97) and equation (109) we obtain, 


ð 18 |/dah\ ə 3 dQ 
dtr dr (2) ðr (ver E) may 


This form is valid for an arbitrary (fixed) profile of an- 
gular velocity and angular momentum in the disk. Very 
often we are interested in the case of a disk that orbits a 
Newtonian point mass M. In that limit, 


Q = Qk = VGM/r3, (111) 
h = VGMr. (112) 
The radial velocity is given by equation (109) as, 
-3 d 1/2 
Ur = -S2 dr (ver ) 5 (113) 
and the evolution equation has the form, 
OX 38 | 479 1/2 
ðt =r Or b Or (vr ) i (a) 


The surface density thus evolves according to a diffusive 
partial differential equation, whose precise character de- 
pends upon the nature of the viscosity. The equation is 
linear if v Æ f(X), though there is no general reason for 
this to be the case. 


2. Steady solutions 


Steady solutions to equation (114) are easily derived. 
Setting 0X /ðt = 0 and integrating we find that, 


vE = c + cort., (115) 


Determining the constants of integration takes a little 
more work, and the right answer depends on the physics 
of the disk being modeled. It’s easiest to start from equa- 
tion (109). Noting that the accretion rate M is, 


M = —2rr£ū,, (116) 
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and that M must be constant for a steady solution, we 
integrate equation (109). The result is, 
M 


/ dQ 
h = vr? — + const. 


27T dr me) 


For Keplerian (point mass) forms for h and Q, the con- 
stant term is negligible at large r. We have, as r — oo, 


vo (118) 


30 
The surface density of the disk is inversely proportional 
to the viscosity. 

To get at the second constant of integration, we note 
that the constant appearing in equation (117) is propor- 
tional to an angular momentum flux Mh. To obtain the 
standard form of the steady disk solution we assume that 
at some radius r = ř the first term on the right-hand-side 
of equation (117), which is proportional to the viscous 
torque G, vanishes. If h and Q are given by Keplerian 
expressions, we then find, 


v% = a i — V" F 
37 r 
This is the solution for a steady-state disk subject to a 
zero-torque boundary condition at r = ř. Classically, 
this boundary condition can be physically justified for a 
disk around a slowly rotating, non-magnetized star, with 
T ~r,, the stellar radius (Pringle, 1977), and for a disk 
around a black hole, in which case f can be identified with 
the innermost stable circular orbit (Bardeen, 1970; Page 
and Thorne, 1974; Shakura and Sunyaev, 1973). How- 
ever, in neither situation is the justification watertight 
(for the black hole case see, e.g.; Agol and Krolik, 2000; 
Gammie, 1999a). 
The heating rate per unit volume in the disk is given 
by, 


(119) 


do dQ\* 9 
= Tror — = —} =>? 120 
ae ror ap u(r) 4” K 5120) 
where the last equality applies only for a point-mass Ke- 
plerian potential. Integrating over z, the heating rate per 
unit surface area in the disk plane is, 


~ 9 9 
Q4 =| {lOicdz = rae (121) 
This heat may result in an increase in the temperature 
of the disk, and it may be transported radially by the 
disk flow. If these effects are negligible and the energy is 
radiated locally, the disk effective temperature is, 
9 

20T ae = qua, (122) 
where ø is the Stefan-Boltzmann constant and the factor 
of two comes from the fact that the disk radiates from 


FIG. 10 Geometry for computing the radial temperature pro- 
file of a disk primarily heated by irradiation from a central 
source, rather than by internal dissipation. 


both its upper and lower surfaces. This equation does 
not require that the disk be in a steady state. If the disk 
is in a steady state, however, with the profile given by 
equation (119), the temperature distribution is, 


3GMM 


Th = 
= 8ror? 


e 


r 


"| . (123) 


The steady state temperature profile does not depend on 
the viscosity. 

The temperature profile for a steady viscous disk is 
not what you get from a fully local toy model in which 
available gravitational potential energy is lost as radi- 
ation at every radius. To see this, suppose that mass 
Am at radius r moves to r — Ar while remaining on a 
near circular orbit. Half of the liberated potential energy 
goes into increased kinetic energy, so the energy avail- 
able to heat up the gas is GMAmAr/2r?. If the time 
scale for mass to move inward is At, and the heat is 
radiated uniformly from the annulus with total surface 
area 4rrAr, the expected temperature profile would be 
T4; = GMM /(8ror?). At large r, this differs by a fac- 
tor of three from the disk profile given by equation (123). 
Viscous torques cause a significant radial redistribution 
of energy in accretion disks. 


3. Irradiated disks 


The temperature profile given by equation (123) ap- 
plies if the dominant source of disk heating is internal 
dissipation. It is also possible for the dominant source 
to be external irradiation, for example if the accreting 
object is a star. The temperature profile in this limit 
depends upon the shape of the disk, which determines 
the fraction of stellar photons that are absorbed at each 
radius. The simplest case is a flat disk in the equato- 
rial plane, that absorbs all the incident stellar photons 
and re-emits the energy locally as a single temperature 
blackbody. 

To compute the resulting Tog(7), consider a surface in 
the plane of the disk at distance r from a star of radius 
R,. The star is assumed to be a sphere of constant bright- 
ness Jx. Setting up spherical polar coordinates such that 
the axis of the coordinate system points to the center of 
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the star, as shown in Figure 10, the stellar flux passing 
through this surface is 


F= le sin 0 cos dQ, (124) 


where dQ represents the element of solid angle. We count 
the flux coming from the top half of the star only (and 
equate that to radiation from only the top surface of the 
disk), so the limits on the integral are, 


—1/2<o< 1/2 


0 << sin! (=) . 
r 


Substituting dQ = sin 0dðd¢, the integral for the flux is, 


(125) 


T/2 sin™t(R«/r) 
F=], / cos do | sin? @d0, (126) 
0 


—n/2 


which evaluates to, 


R, R, R,\? 
F=], sin ( ) ( ) 1 ( ) (127) 
r r r 
For a star with effective temperature Tą, the brightness 
I, = (1/r)oT?, with o the Stefan—Boltzmann constant 
(e.g. Rybicki and Lightman, 1979). Equating F to the 


one-sided disk emission oT4, we obtain a radial temper- 
ature profile, 


Th l i À A 
ef Nig 2 fa (128) 
T+ x r r r 


Integrating over radii, we obtain the total disk luminosity, 


Lis = 2 X | QrroT dr 


bk 


= LAR 
A flat disk that extends all the way to the stellar equator 
intercepts a quarter of the stellar flux. 

The temperature profile given by equation (128) is 
approximately a power-law at large radii. Expanding 
the right-hand-side in a Taylor series in the limit that 
(R../r) <1 (i.e. far from the stellar surface), we obtain, 


(129) 


Tec, (130) 


as the limiting temperature profile of a thin, flat, passive 
disk. For fixed molecular weight u this in turn implies a 
sound speed profile 


cs x 78/8. (131) 


and a scaling of the geometric thickness with radius, 


(132) 


An irradiated disk therefore flares (i.e. has a concave 
shape) toward larger radii. If the disk does flare then the 
outer regions intercept a larger fraction of stellar photons, 
leading to a higher temperature. As a consequence, a 
temperature profile Tog œ r—?/4 is the steepest profile 
we would expect to obtain for a passive disk. 

Irradiation is frequently important for protoplanetary 
disks, with standard models (that include self-consistent 
treatments of disk flaring) having effective temperature 
profiles close to Tog œ r—!/? (Chiang and Goldreich, 
1997; Kenyon and Hartmann, 1987). It can also be im- 
portant in high energy accretion environments, for exam- 
ple in X-ray binaries where irradiation of the outer disk 
by X-rays from the inner disk can dominate the local 
thermal balance (Dubus et al., 1999). 


4. Green's function solution 


Assume for simplicity that the viscosity v(X,7r,...) is 
a constant. The surface density of the disk U(r, t) then 
obeys the equation, 
Od 3v ð 1/2 0 1 
~ m), 133 
ot r Or i rat á 138) 


To solve this equation we first manipulate it into the 
standard form of a Bessel’s equation?. We look for a 
solution in which the variables are separated, and modes 
have a decaying time dependence, 


X(r, t) = r8o(r) exp[—A]. (134) 


Here à > 0 and by writing the spatial dependence as 
r®a(r) we have given ourselves a free parameter in £. 


Substituting, 
3v d 1/2 d ( Pk) 
— b — (or ; 


After evaluating the derivatives and dividing through by 


r?—2 we have, 
d?o 3\ do 1 AÀ 
2 i, = Se ke Spe 
ra + (26+ 3)r2 +5(94+5)o+ er oa =0. 
(136) 


Defining k? = A/(3v) and using the freedom to choose 
p = —1/4, 


d?o do 1 
24 0 ao 22 + z 
Petrit (er 5e 0. 


— àro = (135) 


(137) 


3 We’re working toward the famous solution found by Lynden-Bell 
and Pringle (1974), but here following Gordon Ogilvie’s notes on 
“Accretion Disks” from Part III of the Cambridge Mathematical 
Tripos. The solution strategy is still not all that obvious, though 
you might note that we have a diffusion equation in cylindrical 
co-ordinates, which is analogous to a classical example of Bessel’s 
equation—heat diffusion in a cylinder. 
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This is in the form of Bessel’s equation, which has a gen- 
eral solution, 


o = cı Jı /4(kr) + c2Y1/4(kr), (138) 


where cı and cz are constants and Jıj4 and Yj/4 are 
Bessel functions of the first and second kinds respectively. 
The term involving Y; /4(kr) implies a non-zero torque as 
r — 0, so in the case of a point mass that does not spin 
up the disk material c2 = 0. The solution is therefore, 


ux r—V/4 J, ja(kr) exp[-3vk’t]. (139) 


The properties of Bessel functions allow us to write a 
general initial condition for the surface density in the 
form, 


d(r,0) = f > g(k)r—/4 Jy ja(kr)dk. (140) 


The time-dependent solution will then be, 
Se= | glk)r™ 4J 4(ker) exp[—3vk2t]dk. (141) 
0 


The problem is thus solved provided that we can deter- 
mine the decomposition of the initial surface density into 
Bessel functions, given by g(k). 

To determine g(k) we make use of the Fourier-Bessel 
(or Hankel) transform pair. The textbook definition of 
this pair is, 


y= fo f(r) Im (kr)rdr, (142) 
f(r= [ g(k) Im (kr )kdk. (143) 
Writing equation (140) in this form, 
r/4d(r,0) = i j k~1 g(k) Ji ya(kr)kdk, (144) 
the inverse transform is, 
k-!g(k) = a s"4X(s,0)Jıja(ks)sds. (145) 


Substituting in equation (141) the general solution is, 
p—1/4 [ [> 


We express this in the form, 


D(r,t) = 2 G(r, s,t)X(s,0)ds 


(s, 0) Jaya (ks) Jy /4(kr) 


x exp[—3vk7t]s°/4+kdsdk(146) 


(147) 
where, 


G(r, s,t) = Siet] Jıjalks)Jıja(kr) 
0 


kexp[—3vk?t]dk, (148) 
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FIG. 11 The evolution of a thin ring of gas, initially at r = ro, 
under the action of a constant viscosity (equation 153). The 
curves are plotted at dimensionless times that are multiples 
of two, starting at 7 = 0.008 and going up to T = 0.512. 


is the Green’s function. This integral evaluates to, 


p—1/455/4 rs r2 + 8? 
G(r, s,t) = hya( Je | ( i ; 


6vt 6vt 12vt 
(149) 
with I,/4 being a modified Bessel function. It is illus- 
trative to consider the solution for an initial ring of gas 
orbiting at s = ro. Taking the initial condition as, 
m 


X(s,0) = 


6(8 — ro), (150) 


2779 
with ô being a Dirac delta function, the solution follows 
immediately from equations (147) and (148). It can be 
written compactly in terms of dimensionless variables x 
and 7, 


t=, (151) 
ro 
12 
r= (152) 
ro 
In terms of these variables, 
m «4 2x 14+ 2? 
(z, T) = — ha ( ) ex -12 ; 
nr rT T T 
(153) 


The solution described by this equation is shown in Fig- 
ure 11. It displays asymmetric diffusion, with mass flow- 
ing toward r = 0 while the conserved angular momen- 
tum is carried by a vanishing fraction of the mass toward 
r= œ. Although derived under the restrictive and nor- 
mally unrealistic assumption that v is constant, these 
properties are qualitative features of viscous disk evolu- 
tion in the case where there is zero-torque at the inner 
edge of the disk. 

Although we will not discuss the details here, time- 
dependent solutions can also be derived that dispense 
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with the zero-torque assumption (Nixon and Pringle, 
2020; Rafikov, 2016). Particularly simple decretion disk 
solutions (aspects of which were already discussed in §2.5 
of Lynden-Bell and Pringle, 1974) exist if one assumes 
that an external torque maintains a v, = 0 boundary 
condition at a finite radius rj, (Pringle, 1991). It is also 
possible to derive a relativistic version of the disk evolu- 
tion equation (Balbus, 2017). 


5. Self-similar solution 


Lynden-Bell and Pringle (1974) also derived a self- 
similar solution to the disk evolution equation (equa- 
tion 114), for the case where the viscosity is a power-law 
function of radius, 

vor’. (154) 
If a disk with characteristic size rı at t = 0 has a surface 
density profile of the form, 


Z(t =0) = 


3a fI 


exp [-#-) , (155) 


where C is a constant, * = r/r,, and vı = v(r,), then 
the time-dependent solution is, 


C p2-7) 
NFT) = — T7 5/2-9)/ (2-7) = 56 
@,T) =~ exp | (156) 
where, 
T= + Fi (157) 
1 r? 
ts = = 158 
T a 


This solution has proved to be quite useful for comparing 
theoretical models of viscous disk evolution to data (e.g. 
Hartmann et al., 1998). It can be generalized to the case 
where disk evolution is driven by a combination of viscous 
transport and MHD winds (Tabone et al., 2021). 


B. The a-prescription 


A predictive model for disk evolution follows from 
equation (114) if we can write down how the stress, or 
equivalently the viscosity, depends on properties of the 
disk. Shakura and Sunyaev (1973) advanced physical ar- 


guments in favor of the form, 
Tro = —aP, (159) 


where P is the pressure and a is a dimensionless parame- 
tert. This ansatz is known as the “a prescription”. Using 


4 Shakura and Sunyaev (1973) first introduce a with an expression 
involving the magnetic field (H in their notation), similar but 


equation (106) and writing v = p/p, for a Keplerian disk 
an equivalent form is, 
v= zech ~ acsh. (160) 

In keeping with the approximate nature of this exercise, 
the factor of two-thirds is typically ignored and the a- 
prescription written as just V = ac,h. 

To order of magnitude, the microphysical viscosity of 
a fluid can be written in terms of the thermal velocity of 
the molecules vin and the mean-free-path l as, 


(161) 


V ~ Ue pl. 


By analogy, one can view equation (160) as describing an 
effective viscosity due to turbulent eddies whose speed 
scales with the sound speed and whose size scales with 
the disk thickness. Since turbulent velocities exceeding 
the sound speed would cause shocks and rapid dissipa- 
tion, and isotropic eddies could not significantly exceed 
h, this argument bounds a < 1. This argument is not 
terribly useful, as physical mechanisms for disk turbu- 
lence do not yield turbulent structures that look much 
like eddies on scales of the order of h. It is more instruc- 
tive to follow Shakura and Sunyaev’s original train of 
thought, and express a in terms of fluid (Reynolds) and 
magnetic (Maxwell) stresses in a turbulent fluid (Balbus 
and Hawley, 1998), 


(162) 


where the angle brackets indicate a density weighted av- 
erage over space and time. The first term is the Reynolds 
stress from correlated fluctuations in the radial veloc- 
ity and perturbed azimuthal velocity, the second is the 
Maxwell stress from MHD turbulence. 


C. Time scales 


For a thin disk we can express several relevant time 
scales as simple functions of the disk properties. The 
dynamical time scale, 


1 


tayn = Q (163) 


is evidently 1/27 of the orbital period. The time scale for 
establishing vertical hydrostatic equilibrium is the sound 
crossing time across a scale height, thyaro ~ h/cs Using 
h = cs /Q we have, 

h 1 


thydro ~Y — ~ = ~ tayn- 
y Cs Q y: 


(164) 


not identical to our equation (162). The famous version is their 
equation (1.2), T.g = —apc?, where the sound speed includes 
contributions from both gas and radiation pressure. 
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Vertical hydrostatic equilibrium is thus established in a 
circular disk on the shortest possible time scale. In an 
eccentric disk, however, where the gravitational potential 
experienced by a fluid element varies around the orbit, 
the approximate equality between thydro and tayn means 
that hydrostatic equilibrium is not established and inter- 
esting coupled dynamics between the radial and vertical 
structure is possible (Ogilvie and Barker, 2014). 

The thermal time scale is the time scale on which the 
disk would cool if heating processes were instantaneously 
cut off. The thermal energy per unit surface area of 
the disk is U ~ Nc?. Using equation (121) and equa- 
tion (160), 


U 1 


E Go (165) 


tth 
The thermal time scale is the shortest time scale on 
which we expect the emission from an annulus of the 
disk, heated by viscous-like dissipation, to change. 

The viscous time scale is the time scale on which redis- 
tribution of angular momentum leads to gas inflow. If the 
surface density is not in a steady-state, it is also the time 
scale over which the surface density evolves. Starting 
from equation (114), we can estimate the viscous time 
scale by writing the evolution equation in a form that 
more closely resembles a prototypical one-dimensional 
diffusion equation. Defining, 


X = 21/2, (166) 
J= i (167) 


and assuming that v is a constant, the evolution equation 
is, 


Of af 
BE” Daye (168) 
with diffusion coefficient D given by, 
12v 


The diffusion time scale across scale AX implied by equa- 

tion (168) is (AX)?/D. Returning to the original vari- 

ables, the time scale over which viscosity smooths out 

surface density gradients on radial scale Ar is, 
(Ar)? 


Tvisc p : 


(170) 


If the disk has size r, the surface density can evolve on a 
time scale, 


r2 


Tyisc ~ —. 171 
isc y ( ) 
Using the a-prescription, we obtain, 


"i i /h\? 
V1SC aQ r * 


(172) 
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FIG. 12 Estimates of the dynamical, thermal and viscous 
time scales at r = 10GM/c? for disks around supermassive 
black holes of different masses. The disk aspect ratio is taken 
to be h/r = 0.1. The colored bands for the thermal and 
viscous time scales show the range obtained assuming 107° < 
a<0.1. 


For a thin disk the viscous time scale is substantially 
longer than the thermal time scale, and absent special 
circumstances thermal equilibrium is maintained in the 
vertical direction while the surface density evolves due to 
accretion. 

The time scale hierarchy, 


tdyn ee thydro <K tth X tise, (173) 


is a generic property of geometrically thin disks, and is 
one of the main reasons why thin disk theory is internally 
consistent and useful (e.g. Pringle, 1981). Figure 12 gives 
a concrete example of these time scales at ten gravita- 
tional radii around supermassive black holes of various 
masses, for h/r = 0.1 and different assumed values for 
a. For a 10° Mo black hole, for example, the dynamical 
time scale in the inner disk is of the order of a day, the 
thermal time scale is around a month, and the viscous 
time scale is around ten years. 


D. a-model disks 


Adopting the a-prescription (§III.B) the dependence 
of v on the local disk conditions and on a, v(a, X, Qx), 
can be determined. With this function in hand, equa- 
tion (114) can be solved (usually numerically) for the 
time-dependent evolution of an arbitrary initial surface 
density profile. Steady-state solutions (usually analytic) 
for the surface density profile ©(r,M) can also be found. 
These are “a-model” or “Shakura-Sunyaev” disk solu- 
tions. Recall that none of this effort is necessary if our 
only concern is the profile of the disk effective tempera- 
ture in steady state, as that is given by equation (123) 
independent of the form of the viscosity. 
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A toy example shows how a-model disks are con- 
structed. Assume, for no reason other than simplicity, 
that the vertical structure of the disk is isothermal. The 
effective temperature Tig is then the only temperature 
characterizing the disk at some radius, and the viscos- 
ity can be derived from a triplet of already-introduced 
equations, 


V = QCsh, (174) 
Cs 
h= 1 
A (175) 
9 
20T = que (176) 


The sound speed is related to the temperature through, 
2 _ kpTor 


(177) 


Cc ? 
S UMp 
where kg is Boltzmann’s constant, mp is the mass of the 
proton, and u is the mean molecular weight in units of 
Mp. Using this, we eliminate Tog from equations (174)- 
(176) and obtain an expression for the viscosity, 


1/3 k 4/3 
9 (=) ag PDS, 


= 1 
=e oe (178) 


For fixed central mass M, the predicted viscosity scales as 
v x a4/37y)!/3_ The equation for the evolution of the disk 
surface density, equation (114), would be non-linear with 
this viscosity. In steady-state, vi oc M (equation 119), 
so away from the inner boundary the disk surface density 
would scale as © œx M3/4r—3/4, We have not specified a, 
but as long as this can be taken to be fixed (determined, 
perhaps, from simulations of physical angular momentum 
transport mechanisms or inferred from observations of 
time-dependent disk systems) we have a full solution for 
the evolution of geometrically thin disks. 

The toy model given above captures the spirit of a- 
model disks, but it’s not quite the full story. Real disks 
will not be vertically isothermal. This extra complexity 
can be captured at different levels of approximation, 


(i) At each radius, characterize the disk’s vertical 
structure in terms of a central temperature Te as 
well as an effective temperature Tog. We can derive 
a relation between Te and Teg by considering how 
energy is transported vertically within the disk. 
Then, assuming that the sound speed that enters 
into the expression for the viscosity (v = ac?/OQx) 
corresponds to the central temperature, we can pro- 
ceed as before and derive the functional form of the 
viscosity. This is described as a “one-zone” model 
for the disk vertical structure. 


(ii) Alternatively, we can write down and solve (numer- 
ically) differential equations for the vertical disk 
structure, p(z), T(z), in a manner directly anal- 
ogous to the equations of stellar structure. This 


approach requires a point-by-point specification of 
the stress, for which we could adopt the original 
Shakura and Sunyaev (1973) form, T,4 = —aP, or 
something else perhaps derived from simulations. 
Because of the separation between the thermal and 
viscous time scales in a geometrically thin disk, it is 
normally consistent to solve for the vertical struc- 
ture separately from the radial structure. This is 
sometimes called a “1+1D” disk model. 


Both of these approaches are well-defined. Whether the 
additional complexity of a vertical integration leads to 
a physically more realistic model, however, is an open 
question. Unlike in the case of stellar structure—where 
the energy transport processes and rate of nuclear en- 
ergy generation are quite well-known—the physical pro- 
cesses entering into calculations of disk vertical structure 
are uncertain. I wouldn’t ascribe much physical reality 
to modest differences between disk models with differing 
formal degrees of approximation. 

To write down a version of the one zone equations, suit- 
able for deducing the properties of steady-state a-disks, 
we need only a relation between the central temperature 
Te and the effective temperature Teg. As in stellar struc- 
ture, energy can be transported from the hot interior to 
the cooler photosphere by radiative diffusion or by con- 
vection®. Consider the limit of an optically thick disk 
with radiative transport. The vertical energy flux is (e.g. 
Rybicki and Lightman, 1979), 


lol? dT „4 


F = — —_ — = 


(179) 
where KR is the Rosseland mean opacity (with units of 
cm? g~'). In equating the flux to a constant we have 
assumed that energy dissipation is strongly concentrated 
at z = 0. Noting that the increment of optical depth 
dr = krdz, we integrate from the mid-plane to the pho- 
tosphere, 


16 Teff 


Zph 
Para Tig | dr. (180) 
0 


Te 


If the disk is sufficiently optically thick that T, >> Tog 


then we obtain, 
TA 3 
~T, 
Tef 4 


as the relation between the central and photospheric disk 
conditions. 


(181) 


5 Turbulent transport or transport by waves are also in principle 
possible. In fact, they may well be important, and we ignore 
them here only because they are harder to capture in simple 
analytic formulae. 
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With this expression in hand, we can write down a set 
of equations that determine the steady-state radial struc- 
ture of a-model disks in the one zone approximation. The 
disk is specified by the central mass M, accretion rate M, 
and innermost radius ř, where a zero-torque inner bound- 
ary condition is imposed. The variables to be determined 
are the mid-plane density po, pressure P, temperature Te, 
sound speed cs, surface density 4, scale height h, optical 
depth to the mid-plane 7, Rosseland mean opacity kp, 
and viscosity v. We assume an a-prescription, with a 
a constant, and approximate the opacity as a power-law 
function of the central density and temperature. Col- 
lecting together previous results (equations 10, 8, 5, 123, 
160, and 119), and adding in an equation of state together 
with some straightforward definitions, the set is, 


NES 182) 
po = h 
Cs 
h = 1 
On’ 83) 
P 
on =, 184) 
Po 
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P=- t+ Ti, 185) 
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1 
T= ghRy, 189) 
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Up to a few not-so-important numerical factors, these 
are the standard equations used to determine thin disk 
structure in the Newtonian limit. Additional discussion 
of them can be found in Frank et al. (2002). 

As written, the mid-plane pressure is the sum of a gas 
pressure component and one due to radiation pressure. 
Usually one or the other of these pressure sources is much 
larger than the other, with as a rule of thumb radiation 
pressure dominating in black hole disks close to the inner- 
most stable circular orbit, and gas pressure dominating 
otherwise®. Dropping either gas or radiation pressure, 
we can solve the set of equations for the steady-state 
disk structure, verifying after the fact that we dropped 


6 Under some circumstances, specifically when the disk is threaded 
by a net magnetic field, magnetic pressure, Pg = B?/87, may 
also contribute to vertical support against gravity (Bai and 
Stone, 2013; Mishra et al., 2020; Salvesen et al., 2016; Zhu and 
Stone, 2018). This possibility is not normally considered in clas- 
sical disk models. 


the right one. Away from the inner boundary, the solu- 
tions take the form of power-laws, e.g. X œ rY M® Ma7, 
with power-law indices that depend upon the source of 
pressure and upon the functional form of the opacity. Be- 
cause radius and central mass only enter the equations 
combined in the form of the Keplerian angular velocity, 
the solutions only depend on Ox. 


E. Self-gravitating a-disks 


In most cases, and in particular when the source of an- 
gular momentum transport is the MRI, the a-disk equa- 
tions in no way determine a (though the formalism 
would be inconsistent for a values large enough to in- 
duce supersonic inflow). Self-gravitating disks can be an 
exception. Their stability is a function of the Toomre 
Q = ¢c.Qx/(7GD) (ILE, specializing to a Keplerian 
disk). Suppose, somewhat reasonably, that the angu- 
lar momentum transport rate is a function of the local 
disk conditions and increases rapidly as Q drops below 
some critical value Qo ~ 1. Under these assumptions, the 
combination of self-gravitating transport and local ther- 
mal equilibrium can establish a stable feedback loop that 
maintains Q ~ Qo. If Q > Qo, reduced transport leads 
to reduced heating, lowering cs to re-establish Q = Qo. 
The reverse happens for Q < Qo. Self-gravitating disks 
are then expected to be everywhere marginally stable, 
with Q ~ Qo. This type of model was introduced by 
Paczynski (1978). 

Imposing Q = Qo in addition to the usual set of a- 
disk equations has an important consequence: a is no 
longer a free parameter but rather a specified function of 
the local disk conditions (Gammie, 2001; Levin, 2007). 
The evolution of the disk in this limit is then fully de- 
termined. Rafikov (2015), and references therein, detail 
such “gravito-turbulent” disk models. They are useful 
provided that non-local angular momentum transport 
and mass infall (a complication that often accompanies 
self-gravity in astrophysically relevant settings) can be 
consistently ignored. 


F. The values of a 


Although a great deal of effort has been expended over 
the years in trying to determine “the” value of a, it 
should be clear from the discussion so far that this is 
an illusory quest. Even to the extent that œ provides a 
good parameterization of the strength of accretion disk 
turbulence, its value ought to depend upon the MHD 
properties of the disk, on the strength of self-gravity, and 
so on. The consensus theoretical expectation is that for 
disks that are well-described by ideal MHD, turbulence 
in the absence of net magnetic flux yields (Davis et al., 
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2010; Simon et al., 2012), 


aznrp œ 0.01 — 0.02. (191) 
There remains some uncertainty about how well con- 
verged this computational result is (Ryan et al., 2017). 
In the presence of a net vertical magnetic field B, trans- 
port is stronger, with local simulation results indicating 
that (Bai and Stone, 2013; Hawley et al., 1995; Salvesen 
et al., 2016), 

ayp © p7”, (192) 
where p, is defined as the ratio of the gas pressure to the 
magnetic pressure in the net field at the disk mid-plane, 
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This scaling implies that disks with moderately strong 
vertical fields, 8, < 10?, are strongly turbulent and 
generate mid-plane toroidal fields with magnetic pres- 
sure comparable to the gas pressure. Such magnetically 
dominated or “magnetically elevated” disks have stability 
properties that differ in interesting ways from standard 
Shakura-Sunyaev disks (Begelman and Pringle, 2007). 

For self-gravity, local gravito-turbulent models predict 
that a scales with the local cooling time as (Gammie, 
2001), 
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with an upper limit set by fragmentation at agg ~ 0.1 
(Gammie, 2001; Rice et al., 2005). As with the MRI, pin- 
ning down these numbers precisely from numerical sim- 
ulations is none too easy a task. 

Observationally, œ can be estimated in systems where 
variability exposes the viscous time scale (equation 170), 
with the most important example being dwarf novae, 
whose disks show limit cycle behavior due to thermal in- 
stability (see §VII.B). Dwarf nova outbursts can be well- 
described as time-dependent a-disks (Bath and Pringle, 
1982; Meyer and Meyer-Hofmeister, 1981; Mineshige and 
Osaki, 1983a; Smak, 1984). Under well-ionized con- 
ditions, modeling of these systems suggests a ~ 0.1 
(Hameury, 2020; King et al., 2007). The inferred larger 
value of a, as compared to predictions from simplified 
MRI simulations, may be due to convection in dwarf nova 
disks (Hirose et al., 2014). The strength of turbulence in 
protoplanetary disks can be constrained more directly 
from observations of molecular line broadening (Hughes 
et al., 2011). Such analyses suggest that much lower lev- 
els of turbulence (in some cases only upper limits are 
obtained) occur in very weakly ionized disks (Flaherty 
et al., 2020, 2017). 


G. The Shakura-Sunyaev solution 


Thin disk solutions depend upon whether the main 
source of pressure is gas or radiation, and upon the opac- 
ity under conditions encountered near the disk mid-plane. 
For disks around compact objects (black holes, neutron 
stars, and white dwarfs) two opacity regimes cover most 
conditions of interest. At high temperatures, electron 
scattering dominates. For plasma with a typical astro- 
physical distribution of elements, the opacity is, 


Kos = 0.34 cm? g7}. (195) 


At lower temperatures, free-free opacity, which can be 
approximated using Kramers’ law, applies, 


kg = 6.4 x 10?2pT-7/? cm? got. (196) 
(Here, p and T are understood to be expressed in c.g.s. 
units.) Kramers’ law remains valid down to the temper- 
ature where hydrogen recombines, at T ~ 1—2 x 104 K. 
At lower temperatures, which can be encountered at large 
radii in AGN disks and which are typical of protoplan- 
etary disks, molecules, dust and ice grains provide the 
opacity. 

The properties of Shakura-Sunyaev disk solutions are 
not terribly intuitive, but one important result—the disk 
thickness in the innermost radiation pressure dominated 
region—is readily derived. The vertical flux of momen- 
tum carried by radiation, F,/c, is equal to (using equa- 
tion 123), 
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At high enough temperatures the force that the radia- 
tion exerts per unit mass on the gas is F,a7/c, where 
op = 6.7 x 107?° cm? is the Thomson cross-section ap- 
propriate for electron scattering. Setting this equal to 
the vertical acceleration due to the gravity of the central 
object, which in the Newtonian limit is just Q2z, the 
scale height in the radiation pressure / electron scatter- 
ing dominated regime is, 


30: T 
v= 87 | 
Away from the inner boundary, the scale height (h itself, 
not h/r) is predicted to be constant with radius, with a 
value that is proportional to the accretion rate. 
Without further ado, we quote without derivation the 
Shakura and Sunyaev (1973) disk solution, which takes 


the form of piece-wise power-laws corresponding to three 
regimes, 


(197) 


M. (198) 


a 


e An inner disk, dominated by radiation pressure and 
electron scattering opacity. 
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e A middle disk, dominated by gas pressure and elec- 
tron scattering opacity. 


e An outer disk, dominated by gas pressure and free- 
free opacity. 


The extent of the outer disk is limited by the validity 
of the free-free opacity formula—which as noted already 
fails at low temperature—and / or by the onset of other 
physical processes such as gravitational instability. The 
model is conceptually just the solution of equations (182- 
176), though there is some added subtlety that arises 
from the fact that electron scattering is a true scattering 
process that does not alter the energies of either photons 
or electrons. This means that even at high temperatures, 
sub-dominant absorption opacity plays a critical role in 
thermalizing the emergent radiation. There are also or- 
der unity numerical factors that differ between our equa- 
tion set and that of the original Shakura and Sunyaev 
(1973) paper. 

The Shakura and Sunyaev (1973) solutions can usefully 
be expressed in terms of dimensionless variables for mass, 
accretion rate, and radius, 


M 
= — 1 
m Me’ (199) 
M 
n= 200 
m= 3x 10-8(M/Mo) Mo yr!’ vm) 
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where rg, the Schwarzschild radius, is given by, 
2GM 
rs = Z. (202) 
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These scaling are evidently intended for black hole ac- 
cretion problems. The Schwarzschild radius is the radius 
of the event horizon for a non-rotating black hole, which 
has an innermost stable orbit at r = 3rg. The accre- 
tion rate scaling corresponds roughly to the Eddington 
limit, which is also most directly relevant to black hole 
and other energetic accretion environments (see IV.A). 
Nonetheless, these are fundamentally Newtonian solu- 
tions, which can be rewritten with scalings more appro- 
priate to, e.g., white dwarfs, without difficulty. Denoting 
the inner disk with a subscript i, the scale height, surface 
density, and central temperature are, 


hi = 3.2 x 10°rmf cm, 
X; = 4.6a lm ir? fT] g em=?, 


Ta = 2.3 x 107a im 4r? K, 


(203) 
In these expressions, f = (1 — ge) For the middle 
disk, denoted with a subscript m, 

hm = 1.2 x L047 M10 p55 99/10 p!71/?9 1/5 om 

Em = 1.7 x 105074 n3 5 m5r 93/8 g omo?, 


Ta = 3.1 x 10807 1/5 rn?/5m 15r? f2/5 K, (204) 
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FIG. 13 Transition radii rip; (red: between the inner radia- 
tion pressure dominated disk and the middle gas pressure / 
electron scattering dominated disk) and rém (blue: between 
the middle gas pressure / electron scattering dominated disk 
and the outer gas pressure / free-free opacity dominated disk), 
plotted for different black hole masses in the Shakura-Sunyaev 
solution. The solid lines assumes mn = 1, the dashed lines 
m = 0.1, in both cases for a = 0.1. 


For the outer disk, denoted with a subscript o, 


ho =6.1x 103a 1/10 753/20477,9/10,,19/8 3/20 cm, 
Eo =6.1x 105074/5rn7/10m1/57173/4 7/10 g cm~?, 


Teo = 8.6 x 107a 15m3 m15 py! 3/4 3/19 K, (205) 


The transition radii between the regimes (7/,,, for the 


mi 


inner to middle disk, r/m for the middle to outer disk) 


om 
are given implicitly by solving, 
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The dependence on a is weak, so these radii are mostly 
dependent on the accretion rate and black hole mass 
(bearing in mind that the latter can vary across many 
orders of magnitude). 

Figure 13 shows the dependence of r/,, and rim as a 
function of black hole mass. For a stellar mass black 
hole with M = 10 Mo, and m = 1, rl; œ 140 (i.e. 
840 GM/c?) and r/m © 6250. For a supermassive mass 
black hole with M = 10? Mo, again assuming rh = 1, 
ri, 540 and r/m œ 6250. In dimensional units, in the 
supermassive case the transition from radiation to gas 
pressure occurs at about 5 x 10!° cm, while the transition 
from electron scattering to free-free opacity is at about 
6 x 1016 cm (0.02 pe). A lower accretion rate of m = 0.1 
moves both radii in by a factor of 5-6. Qualitatively we 
observe that (1) all three Shakura-Sunyaev regions are 
predicted to be present for reasonable parameter choices, 
and that (2) radiation pressure is relatively more impor- 
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FIG. 14 The geometric thickness of Shakura-Sunyaev disks as 
a function of dimensionless disk radius r’ = r/3rg is shown for 
supermassive (m = 10’, red) and stellar mass (m = 10, green) 
black holes. The solid lines show high accretion rate solutions 
(rn = 1), the dashed lines a lower accretion rate (m = 0.1). 
In all cases a = 0.1. The geometric thickness of the inner 
radiation pressure dominated region depends strongly on the 
accretion rate, and can be quite large for high accretion rates. 
The middle and outer disk regions are much thinner, espe- 
cially for disks around supermassive black holes, where values 
of h/r ~ few x 107° are characteristic. (The small disconti- 
nuities in the plotted curves are due to approximations and 
do not have any physical significance.) 


tant for disks around supermassive black holes as com- 
pared to stellar mass examples. 

Figure 14 shows the radial dependence of the predicted 
geometric thickness of some selected Shakura-Sunyaev 
disk solutions. Although these are “thin” disk solutions, 
the region where radiation pressure dominates (and h is 
constant) is not actually thin at all for m ~ 1. Val- 
ues of h/r significantly in excess of 0.1 are predicted at 
r ~ 10-20 GM/c?. The resulting violation of the as- 
sumptions underlying thin disk theory is remedied in slim 
accretion disk models (Abramowicz et al., 1988), which 
should really be used to give a consistent treatment of 
this region. Conversely, the gas pressure dominated mid- 
dle and outer regions of the Shakura-Sunyaev solution are 
quite thin, especially in the case of supermassive black 
holes. For M > 107 Mo we expect 1073 < h/r < 107? 
in these regions. As a consequence, disk self-gravity 
(SIE, I.E) is predicted to become important for disk 
masses that are much smaller than the mass of the black 
hole. The onset of self-gravity and the likelihood of ensu- 
ing fragmentation, in turn, has far-reaching consequences 
for the radial extent of AGN disks, for the formation of 
stars and compact objects within the accretion disk, and 
for how supermassive black holes are fuelled and grow 
(Goodman, 2003; King and Pringle, 2006; Levin, 2007; 
Shlosman et al., 1990). 

Novikov and Thorne (1973) generalized the Newtonian 
Shakura-Sunyaev thin disk solution to include relativ- 


ity. The Novikov-Thorne solution does not introduce 
any novelties in its treatment of the disk physics, but 
the proper inclusion of all the relativistic effects is just 
as tricky as you might expect. There are plenty of oppor- 
tunities for making mistakes (even the original authors, 
no slouches when it comes to relativity, didn’t get it quite 
right). I recommend Abramowicz and Fragile (2013) as 
a source for the explicit solution, and as a starting point 
for reading the literature on Novikov-Thorne disks. 


IV. ENERGETICS OF DISK ACCRETION 


The thin disk solutions are predicated on two assump- 
tions: that the energy released by accretion can be radi- 
ated on a time scale that is short compared to the local 
inflow time scale, and that the outgoing radiation has a 
negligible impact on the flow dynamics. We now turn 
to what happens when these assumptions fail. Radiation 
becomes dynamically important when the luminosity of 
the disk (or that of the central object) reaches the Fd- 
dington limit. Radiative cooling ceases to be efficient 
both when the accretion rate is very low, due to plasma 
physics effects (the regime of radiatively inefficient accre- 
tion), and when it is very high, due to photon trapping 
(the regime of hyperaccretion). 


A. Eddington limit 


The Eddington limit is the luminosity at which radia- 
tion pressure from a central point source balances gravity, 
curtailing spherically symmetric accretion. Noting that 
photons of energy E carry momentum p = E/c, the mo- 
mentum flux at distance r from an isotropic source with 
luminosity L is L/(4mcr?). If the opacity of the gas is x, 
the outward radiative force per unit mass of gas is, 


KL 
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Equating to the inward force per unit mass of gravity, 


GM 
foray =~) (209) 
r 
the Eddington limiting luminosity is, 
4ncGM 
Lyaa = ————. (210) 


At sufficiently high temperatures the opacity is due to 
Thomson (electron) scattering, and k = or/my, where 
or = 6.7 x 10~?° cm? is the Thomson cross-section and 
my = 1.66 x 10-74 g is the mass of a hydrogen atom. 
(The radiative force acts on the electrons, but these are 
tightly coupled to the protons electrostatically.) Under 
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these conditions, numerically, 


4 GM 
ie. (211) 
OT 

M 
~ 1.2 x 10° (=) erg s7}, (212) 

M 
~ 3.2 x 104| —} Lo. 213 
3.2 x 10 (=) (213) 


A modest correction is needed for Thomson scattering 
on gas that is not pure hydrogen, with a larger one being 
in order if the surrounding gas is cool and the opacity is 
due to dust. Dust opacity can be relevant to accretion 
during massive star formation, and to the dynamics of 
gas at relatively large distances from supermassive black 
holes in AGN. 

It is hopefully obvious that the Eddington limit, de- 
rived assuming spherical symmetry, is best-regarded as a 
characteristic luminosity above which radiative forces are 
guaranteed to matter for surrounding accretion flows. As 
with the Pirate’s Code in Pirates of the Caribbean, the 
Eddington limit is more what you call a guideline than 
an actual rule. 


B. Radiative efficiency of thin disks 


The Newtonian estimate for the luminosity of accretion 
at rate M onto a star with mass M and radius R, is just, 
GMM 


L= : 
Ry. 


(214) 


For a black hole the presence of an event horizon means 
that a fraction of the potential energy liberated by ac- 
cretion can, in principle, be “lost” across the horizon, 
augmenting the mass of the hole as measured at large dis- 
tance. The luminosity associated with black hole accre- 
tion then depends on the mode of accretion, and on some 
uncertain MHD physics close to the hole. The standard 
estimate is based upon the assumption that the black 
hole accretes from a disk that is geometrically thin, with 
a zero-torque boundary condition at an innermost radius 
that is close to the hole. 

The fiducial estimate for the radiative efficiency of 
thin disks relies on various properties of Kerr black holes 
(Kerr, 1963). For a rotating, uncharged, black hole, we 
define the dimensionless spin parameter a, in terms of 
the mass M and angular momentum J via, 


(215) 


The spin is limited to 0 < a, < 1 (or —1 < a, < 1 using 
negative values to denote orbits that are counter-rotating 
with respect to the spin). Amongst the many important 
properties of the Kerr metric are the radius of the event 
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FIG. 15 The horizon radius, and that of the innermost stable 
circular orbit for prograde (solid) and retrograde (dashed) 
orbits, is plotted as a function of the dimensionless black hole 
spin parameter ax. 


horizon, 


ra = (14 via) Z, 


and the radii of the innermost stable circular orbits, 
which differ depending on whether the orbits are co- 
rotating or counter-rotating in the equatorial plane of 
the hole (Bardeen et al., 1972; Shapiro and Teukolsky, 
1983), 


(216) 


GM 
TISCO = [3 + Zo F V6 Z1)(3 + 41 4 22Z2)| >? 
Aig) [a +a)” rs a.) 
Zə = 4/302 + 2?. (217) 


These quantities are plotted in Figure 15. The innermost 
stable circular orbit lies at 6GM/c? for a non-rotating 
(Schwarzschild) black hole, and tends towards GM/c? 
and 9GM/c? for co-rotating and counter-rotating orbits 
respectively as a, — 1. The variation of rrsco with a, 
provides the physical underpinning for various black hole 
spin estimates derived from electromagnetic observables 
(Reynolds, 2021). 

Armed with these results, we can define the fiducial es- 
timate of the radiative efficiency properly. Assume that 
a thin disk extends from large radius, where the mass ac- 
cretion rate is M , down to the innermost stable circular 
orbit, where a zero-torque boundary condition is applied. 
Gas interior to the disk—between risco and ry—is as- 
sumed to neither radiate nor to exert any feedback effects 
on the disk flow. The radiative efficiency 7, defined via, 

L=nMe, (218) 
is then fully specified by the value of the binding energy 
of particle orbits at rigco. Defining the ancillary variable 
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FIG. 16 The nominal radiative efficiency of thin disk accre- 
tion is plotted as a function of the dimensionless black hole 
spin parameter ax. Negative values of a. denote accretion 
from equatorial counter-rotating orbits, positive values ac- 


cretion from equatorial co-rotating orbits. For a. = 0 the 
efficiency n = 1 — \/8/9 = 0.0572, for aą = —1 we have 


n = 1 — 1/25/27 = 0.0377, while for a, = 1 (this value of the 
spin is not quite physically realizable) 7 = 1 — ,/1/3 = 0.423. 


E = 1 — ņ the efficiency is found by solving the equation 
(Shapiro and Teukolsky, 1983), 
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The resulting function 7(a,) is shown in Figure 16. For 
a, = 0 the efficiency is 5.7% (i.e. that fraction of the 
rest mass of the accreting gas is radiated from the disk), 
while as a, — 1 the limiting efficiency is 42%." 


1. Salpeter time 


The results for the Eddington limit and the radiative 
efficiency of disk accretion can be combined to give an 
estimate of how rapidly black holes can grow from thin 
disk accretion. To do so, assume that the accretion rate 
is limited to the value that would yield an Eddington- 
limited luminosity. For a black hole of mass Mpu, 

Mc? e 4nrcmyG 
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Subtracting the rest-mass equivalent of the energy lost 
via radiation, accretion at rate M increases the mass of 


7 Although mechanically it is possible to spin up a Kerr black hole 
arbitrarily close to ax = 1, a black hole spun up by disk accretion 
necessarily consumes disk photons on orbits that counteract the 
spin up. Thorne (1974) found that this effect limits the maximum 
spin to ax ~ 0.998, with a radiative efficiency 7 ~ 0.3. 


the black hole according to, 


Mpuy = (1-7) M. (221) 
A black hole that always grows as fast as it can—at the 
Eddington limit—then obeys, 
dMpu An (1 = n) myG 
= Maa. 
dt n COT 


(222) 


At fixed radiative efficiency (i.e. at fixed spin) the result 
is exponential growth, 


t 
Mguu = Mo exp B ; (223) 
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where Mo is some initial mass and tg, the Salpeter time 
(Salpeter, 1964), is a characteristic time scale, 


1 COT 
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ts = (224) 
For 7 = 0.1, the Salpeter time is about 50 Myr. 

The Salpeter time is a useful characteristic time scale 
for how fast black holes accreting from geometrically thin 
disks grow. You will often see it referenced as part of 
an argument about supermassive black hole formation, 
which goes as follows. We observe luminous quasars, 
likely hosting black holes with masses Mgy ~ 10° Mo, 
at redshifts exceeding z = 7 (Mortlock et al., 2011; Wang 
et al., 2021). If we assume (say) that the seed for such 
a quasar formed at z = 20, the time available between 
z = 20 and z = 7 for it to grow is only approximately 
5.9 x 108 yr. This is 11.6 Salpeter times (for 7 = 0.1), 
while 18 e-foldings are needed to grow from a 10 Mo 
stellar mass black hole to a 10? Me supermassive black 
hole. The observed early growth of massive black holes 
in the Universe appears to pose problems, or at least to 
provide strong constraints on the masses of black hole 
seeds (which may not be stellar mass black holes at all; 
Begelman et al., 2006). 

This argument, in my opinion, is rather tired from 
overuse, and in its strong form requires undue faith in 
the Eddington limit being a strict limit. However, it 
does justify the weaker statement that the existence of 
high redshift quasars requires a high duty cycle of prior 
accretion, at a rate high enough to imply that radiative 
forces are important. 

A separate argument, credited to Soltan (1982), takes 
off from the observation that the mass accumulated in 
supermassive black holes, and the total amount of en- 
ergy radiated from accretion during their growth, are two 
sides of the same coin. If supermassive black holes grow 
primarily via thin disk accretion, the total mass in black 
holes per comoving Mpc? at z = 0 is related to the in- 
tegral of the AGN luminosity function over redshift, and 
an appropriate comparison between the two constrains 
the radiative efficiency (Yu and Tremaine, 2002). 
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FIG. 17 Setup for computing the transverse momentum im- 
parted by an electron-proton encounter in the impulse ap- 
proximation. 


C. Electron-ion coupling in low density plasmas 


The validity of thin disk solutions for black hole ac- 
cretion is bounded above at m ~ 1, both by the onset 
of radiatively driven outflows as the Eddington limit is 
exceeded and by the importance of radial advection of en- 
ergy as h/r becomes larger. In the other direction, thin 
disk solutions remain internally consistent for m < 1, 
but they are not unique. A hot, geometrically thick disk 
solution also exists below a critical accretion rate. The 
physical origin of the hot solution is tied to an asymme- 
try in the microphysics of electron-ion plasmas. Heating, 
due to small-scale dissipation of turbulent energy, gives 
energy predominantly to the ions, while cooling, due to 
processes such as free-free emission and synchrotron radi- 
ation, is much more efficient for electrons. At low density, 
the time scale for Coulomb collisions to transfer energy 
from the ions to the electrons becomes long—in some 
cases longer than the time scale for the gas to accrete on 
to the black hole. Unable to cool efficiently, the accretion 
flow becomes (or remains) geometrically thick (Ichimaru, 
1977; Rees et al., 1982; Shapiro et al., 1976). It resembles 
a torus or a doughnut more than a disk. 


The time scale for electron-proton thermal equilibra- 
tion can be computed to better accuracy than we need via 
elementary methods. Step one is to calculate the trans- 
verse momentum that is imparted when an electron, with 
mass Me, flies by a proton, with mass mp, at velocity ve 
and impact parameter b. The setup is shown in Figure 17. 
We work in the impulse approximation, and compute the 
momentum change along the unperturbed straight-line 
trajectory of the electron. The force perpendicular to 
the direction of motion is, 
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The total momentum imparted in the perpendicular di- 
rection as a result of the encounter is just the integral of 


the force over time, 
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A quick trignometric substitution, or an even quicker trip 
to Wolfram Alpha, and we have the answer, 
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Ap. = (227) 
With this result in hand we proceed to step two, which is 
to sum up the effect of many encounters occurring with a 
range of impact parameters bmin < b < bmax. The mean- 
square momentum change is the lowest order quantity 
that does not vanish. If the number density of particles 


is n, 
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min 


The integral gives us In(bmax/bmin), which goes by the 
name of the Coulomb logarithm and is written as ln A. 
The same quantity crops up, for the same reasons, in 
discussions of relaxation in gravitating point-mass sys- 
tems such as globular clusters. In a plasma, bmax should 
be set to be the Debye length (e.g. Thorne and Blandford, 
2017), 


(229) 


because this is the spatial scale on which charges are ef- 
fectively screened. The lower limit, bmin, is either the im- 
pact parameter where a single deflection A0 ~ 1 radian 
(as in the gravitational case), or the particle’s de Broglie 
wavelength, whichever is the larger (Spitzer, 1962). All 
this can be calculated properly as a function of density 
and temperature, but for our purposes the details are 
neither very interesting nor very important. The answer 
depends on bmax/bmin only logarithmically, and taking 
ln A ~ 20 is a good enough approximation. 

For the final step, we switch from thinking about mo- 
mentum to thinking about energy. The electron, with 
kinetic energy (1/2)mev?, changes the energy of the pro- 
ton by an amount (Ap1)?/2mp in a single encounter. 
Considering all encounters, the equilibration time is, 
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Let us assume that the typical electron velocity is, 
3kpgT,\ te 
i= (==) can 
Me 
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Then, noting that the Thomson cross-section (Rybicki 
and Lightman, 1979), 


8m et 
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and defining a parameter expressing how relativistic the 
electrons are, 


kpTe 
he = ; 2 
=a (233) 
equation (230) becomes, 
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The rather sloppy methods we have employed here get 
us surprisingly close to the answer obtained by more dili- 
gent investigators. A two-temperature electron-proton 
plasma, where the two species have Maxwellian velocity 
distributions, comes to a common temperature on a time 
scale’ , 


V 2T m 
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where 6, = kpT,/mpc?. When be > 0p this differs from 
our rough and ready result by only a modest numerical 
pre-factor. 


D. Radiatively inefficient accretion models 


A geometrically thick disk solution is possible if the 
density is low enough that the time scale for electron-ion 
thermal equilibration exceeds the time scale for gas to 
flow into the black hole. This can be written as an up- 
per limit on the accretion rate, in units of the accretion 
rate that would generate an Eddington limiting luminos- 
ity, that depends only on a. To obtain this result, we 
consider a quasi-spherical but rotationally supported ac- 
cretion flow, such that n = n(r) only but v, is set by the 
local rate of angular momentum transport. The continu- 
ity equation is, 

M = —A4nr?urp. (236) 
Substituting the order of magnitude estimates, vp = 
—r/tyise and tyise ~ 1/(@QK) (remember that h/r ~ 1, 
by assumption), and using p = nmy, gives, 


M 
4nrmyoaNkr? 


n(r) (237) 


8 In the plasma context this result is usually cited as Spitzer 
(1962), but Spitzer’s textbook passes the buck to his paper on 
stellar dynamics (Spitzer, 1940) for the details of the calcula- 
tion. That’s fine because, up to the question of the correct value 
of In A, the gravitational and plasma calculations are the same. 


We now substitute this expression for n into equa- 
tion (235) and impose the physical condition that, 


(238) 


tep 2 tyiso: 


Although the electrons cool to lower temperatures than 
the protons, they’re also much less massive, so we can 
assume that 9. > 0, and drop the latter. The inequality 
gives, 


Fa 
M< vV2r 4rmyGMa* (mp 93/2, (239) 
2inA OTC Me 


which is approximately as clear as mud. The expression 
simplifies when we note that many of the quantities on 
the right-hand-side also appear in the formula for the Ed- 
dington limit (equation 213). Let us define an Eddington 
mass accretion rate via, 


Lgaa = NMeaae’, (240) 
such that, 
4 GM 
g= EE, (241) 
COT 


(Note that some authors define Mgaa without the factor 
of 7, or fixing 7 = 0.1.) The inequality then simplifies to, 


M < |5 (Mp 93/22 
Mea YV 2InA\m/ $ : 


It’s not easy to give an elementary estimate for 0e, but 
adopting reasonable values for the parameters, 0. = 0.1, 
7 = 0.1, and ln A = 20, the result is, 


(242) 


M 


; < 0.407. 
Mraa 


(243) 


Our argument has been a simplified version of that given 
in Mahadevan and Quataert (1997), but there are multi- 
ple routes to this result. If a ~ 0.1, the conclusion is that 
a geometrically thick, two-temperature accretion flow, is 
a consistent possibility provided that the accretion rate, 
measured in Eddington units, is below 1073-1072. 

The way in which plasma microphysics and flow macro- 
physics combine to give a simple result for the critical ac- 
cretion rate is quite appealing. It should be obvious that 
our derivation is only good to an order of magnitude, 
but there are also other caveats. The fact that a hot 
two-temperature flow can exist below a critical accretion 
rate does not mean that it must—in principle a denser 
thin disk could be a valid solution below the critical rate. 
There has also been sporadic discussion over many years 
as to whether the plasma physics fundamentals underpin- 
ning the argument—that dissipation primarily heats the 
ions (Gruzinov, 1998; Quataert, 1998), and that Coulomb 
coupling is the sole channel for ion-electron energy trans- 
fer (Begelman and Chiueh, 1988)—are secure. The re- 
view by Yuan and Narayan (2014) is a good place to 
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start for a discussion of these questions, and the extent to 
which the uncertain answers impact astrophysically inter- 
esting conclusions about hot accretion flows. A promising 
development is that the basic plasma physics involved is 
increasingly amenable to direct simulation using particle- 
in-cell techniques (Schekochihin et al., 2019; Sironi and 
Narayan, 2015; Zhdankin et al., 2019). 


1. The limiting ADAF solution 


The microphysical inability of a low-density plasma to 
cool efficiently has two important consequences for the 
structure of low m accretion flows. The plasma, or at 
least the ions within it, becomes hot, and hence geomet- 
rically thick. The larger value of h/r, in turn, means 
that the inflow velocity vy becomes larger than it would 
be in a thin disk. These two properties mean that radial 
advection of heat is important, and can be dominant, for 
the structure of low accretion rate accretion flows. The 
moniker Advection Dominated Accretion Flow (ADAF) is 
sometimes used as generic term for radiatively inefficient 
accretion flows, though it can refer to the specific disk 
model introduced by Narayan and Yi (1994, 1995a,b), 
whose work highlighted the key role of advection and 
kicked off a resurgence of interest in the properties of low 
accretion rate disks. 

In general, geometrically thick disks require at least a 
two-dimensional spatial description (axisymmetric seems 
like it ought to be a reasonable simplification, with flow 
variables being functions of r and polar angle 0), and 
of course they may be time-dependent. Radial pressure 
gradients can not be ignored, the angular velocity need 
not be Keplerian, and the energy equation must include 
the advective term. This all adds up to many additional 
degrees of freedom than we have for thin disks. We may 
abandon all hope and turn to numerical simulations, but 
before reaching that point we can make some progress via 
aggressive simplification of the problem. Narayan and 
Yi (1994), in their original paper on ADAFs, assumed 
that the flow was axisymmetric, steady, and could be de- 
scribed via height-integrated equations rather than a full 
(r,0) treatment (height-integrated here implying some- 
thing more akin to averaging over spherical rather than 
cylindrical shells). The governing equations then read, 


d 
Jr (prhv,) = 0, 244) 
du, 2. 2 1d 2 
Vry Q r = -Nkr rer (pc?) , 245) 
d(Qr?) 1 d 3, dQ 
aa a oh vpr hoy ; 246) 
ds 4 = 
2phu,T — = Q7- Q. 247) 
dr 


The symbols here have the same meanings as in the thin 
disk equations, but we’ve added an energy equation in- 


volving the temperature T, entropy s, and the rates of 
viscous dissipation Q4 and radiative cooling Q7 per unit 
area. Narayan and Yi (1994) go on to derive a self-similar 
solution to these equations for a general y, in terms of a 
parameter f that measures the importance of the advec- 
tive term relative to radiative cooling. 

The essence of the ADAF solution is already present in 
the limiting case where there is no radiative cooling what- 
soever. This limit is amenable to a rather simple analysis. 
Following Blandford and Begelman (1999), we consider a 
quasi-spherical accretion flow with angular velocity Q(r), 
density p(r), and isothermal sound speed c,(r). Angular 
momentum transport is assumed to be inefficient enough 
that vp < Qr. In a steady-state flow, with mass accre- 
tion rate M and inward-directed angular momentum flux 
F, we have, 


—Arr7u,p = M, 
Mr?Q-G = Fi. 


(248) 
(249) 


As before, G is the torque that the disk at radius r exerts 
on the disk exterior to that point. Whereas for a thin 
disk we imagine thin annuli interacting through an effec- 
tive viscous process, here we imagine coupling between 
spherical shells. 

Another conserved quantity is derived from the energy 
equation. The total energy density has contributions 
from potential energy (assumed Newtonian), kinetic en- 
ergy (by assumption of slow inflow the relevant velocity 
is Qr), and enthalpy. The isothermal sound speed cs is 
given in terms of the enthalpy h by, 

a (x-1) h 


S ; 
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(250) 


where y is the effective ratio of specific heats. If non- 
relativistic ions dominate the energy density, we’d expect 
y ~ 5/3 as usual for an ideal gas. Radiation pressure 
would give y = 4/3, though this is not a limit that is so 
relevant to very low m flows. Accretion of gas containing 
energetically dominant small-scale magnetic fields would 
also give y ~ 4/3. Leaving y as a parameter for now, the 
outward-directed energy flux Fg is, 

GM y 


1 
Q 2Q? 
g (Sr r Tip al 


2) M=Fpg. (251) 


This expression has a viscous term and an advective term, 
but omits cooling which is assumed to be neglibly small 
by comparison. 

We now assume that the accretion flow is self-similar, 
which means that it has no preferred scale and “looks 
the same” when rescaled by a multiplicative factor. A 
black hole accretion flow has an inner boundary that is 
set by the radius of the event horizon / innermost stable 
circular orbit, and it must be fed at an outer scale by 
some physical process. These realities necessarily break 
self-similarity. If the inner and outer radii are separated 


30 


by many orders of magnitude, however, it is reasonable to 
think that a flow that cannot cool at all ought to be well- 
approximated as self-similar across much of that radial 
interval’. The self-similar assumptions implies that the 
fluid variables scale as, 


Q(r) x 79/2, (252) 
Ur(r) œ cs(r) œx r712. (253) 
3/2 


The continuity equation then implies that p(r) x r7 
and that the pressure P = pc? « r—5/?. This allows us to 
simplify the radial momentum equation (245). Dropping 
the v,(du,/dr) term because vy < Qr, and inserting the 
self-similar scaling for the pressure, 


oq? GM 52 


(254) 
With these simplifications the equations left to work with 
are (249), (249), (251) and (254). 

We next observe that the term Mr2Q on the left-hand 
side of equation (249) scales as r!/? (as must G), and 
that similarly the terms on the left-hand side of equa- 
tion (251) scale as r~!. These equations can only be sat- 
isfied for arbitrary ranges of r if the right-hand sides van- 
ish, F; = Fg = 0. Eliminating G between equation (249) 
and equation (251) we have, 


1 M 
r2Q? 4 G Y 


2 
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= 2 
Tai C (255) 
Using equation (254) the solution for the angular velocity 
and sound speed is, 

ee 6(y —1) GM 

s (9-5) r’ 
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(97-5) © 
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The difference between this solution and the thin disk 
solutions we have discussed previously is clear. The an- 
gular velocity of ADAF solutions can be strongly sub- 
Keplerian, indeed Q > 0 as y > 5/3. 

To complete the solution G has to be specified. This 
can be done in various ways that are all consistent 
with the spirit of the Shakura-Sunyaev a-prescription. 
Blandford and Begelman (1999), for example, take G = 
4rr?aP, which can be seen as defining a for this type 
of flow. Other choices are possible too (Narayan and Yi, 
1994). Once G is fixed, the density and radial velocity of 
the accretion flow are readily derived. The solution gen- 
eralizes to the case where radiative cooling is a fixed but 


9 The same is obviously not true at all for thin disks. As you can 
see from Figure 14, radial variations in the dominant sources of 
pressure and opacity grossly break self-similarity. 


non-zero fraction of the viscous dissipation rate (Narayan 
and Yi, 1994)?°. 


2. Variations on a theme of radiatively inefficient accretion 


The governing equations from which the ADAF so- 
lution is derived hard-wire assumptions about the ra- 
dial constancy of the mass flux, and the nature of radial 
transport of angular momentum and energy, that can be 
questioned. It is therefore possible to construct alternate 
models for radiatively inefficient accretion flows, that are 
based on different physical principles. One such alter- 
native is motivated by the observation that the ADAF 
solution is unstable to radial convection (see for example 
the simulation in Figure 18). Convection in thick disks 
transports energy outward, but angular momentum in- 
ward, opposing viscous transport. If convection is dom- 
inant, these properties suggest that the governing prin- 
ciple that determines the flow structure is marginal sta- 
bility against convective instability (Narayan et al., 2000; 
Quataert and Gruzinov, 2000). A Convection Dominated 
Accretion Flow (CDAF) has a radial density profile that 
scales as p(r) x r~'/?, substantially shallower than the 
ADAF’s p(r) x r73/2. A second alternative is moti- 
vated by the fact that the supposedly accreting gas in the 
ADAF solution is so hot that it is, at best, weakly bound 
to the central object. The Adiabatic Inflow-Outflow So- 
lution (ADIOS) model assumes that the disk responds 
by driving an outflow from its upper and lower surfaces 
at all radii, that co-exists with inflow near the equator 
(Blandford and Begelman, 1999). M is then no longer a 
constant, but an increasing function of radius. In ADIOS 
models very little of the gas that is supplied to the system 
at large radii ever reaches the vicinity of the black hole. 
A third possibility, though one more typically studied 
with simulations rather than self-similar models, is that 
the magnetic field is strong enough to directly dictate the 
flow dynamics, at least close to the black hole. This is 
the regime of Magnetically Arrested Disks (MAD) (Igu- 
menshchev, 2008; Tchekhovskoy et al., 2011). Analytic 
arguments suggest that geometrically thick disks trans- 
port magnetic flux inward more easily than thin disks 
(Lubow et al., 1994), though the exact conditions that 
lead to MAD solutions and their governing physical prin- 
ciples remain a topic of investigation (Begelman et al., 
2022). 

Which of these physical considerations is most determi- 
native of the structure of radiatively inefficient accretion 


10 Essentially the same equations admit a time-dependent self- 
similar solution for an accretion flow of finite extent (Ogilvie, 
1999b). Although this solution is not often discussed, it’s a use- 
ful way to think about obviously time-dependent problems such 
as accretion that ensues following tidal disruption events. 
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flows in nature cannot be divined by pure thought. Nu- 
merical simulations, and observations of hot flows in sys- 
tems such as the Galactic Center (Marrone et al., 2007) 
and M87 (Event Horizon Telescope Collaboration et al., 
2019), are essential. Yuan and Narayan (2014) and Davis 
and Tchekhovskoy (2020) are good starting points for re- 
views of the current state of the art. 


E. Hyperaccretion 


As we have seen, geometrically thick accretion flows 
arise at low accretion rates because plasma physics effects 
prevent low density gas from cooling efficiently. Thick 
disks can also occur at very high accretion rates, when 
the gas is so dense that photons cannot diffuse away faster 
than they are advected into the black hole. This is the 
regime of hyperaccretion. 

The new aspect of hyperaccreting systems is the exis- 
tence of a trapping radius rtrap, interior to which photons 
cannot escape being dragged into the black hole. Begel- 
man (1978) derived the trapping radius by considering a 
generalized version of spherical Bondi accretion (Bondi, 
1952), which we will get to in §VIII.A. We can obtain the 
main result, however, without needing to specify many 
details of the accretion flow. 

The trapping radius is defined as the location where 
the local infall time matches the time for outward photon 
diffusion. For gas at radius r, with inflow speed vr, the 
infall time is just, 

(258) 


tinfall ~ —. 
To find the diffusion time, we use basic results from 
the theory of random walks (e.g. Rybicki and Lightman, 
1979). If the photon mean-free-path is l, then after N 
scatterings the photon will have diffused an expected dis- 
tance, 


d~ VNI, (259) 
after a time, 
Nl 
Atm ~—. (260) 


The diffusion time scale from radius r is thus, to order of 
magnitude, 
z 
tdiffuse ma Te (261) 
c 
Suppose now that the temperature of the accreting gas 
is high enough that electron scattering provides the main 
source of opacity. The mean-free-path is then, 


jmn (262) 
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FIG. 18 Visualization of a radiatively inefficient accretion flow in the viscous hydrodynamic limit. Such flows are geometrically 
thick, and develop large-scale flow structures that resemble convection. Outflows are likely. (Unpublished simulations; Armitage 


and Dullemond, 2000.) 


where p is the density, my is the mass of a hydrogen 
atom, and or is the Thomson cross-section as usual. The 
density is related to the accretion rate through, 


M =4nr? upp. (263) 


Combining equations, v, cancels out and we find that 
tintan equals taiftuse at a trapping radius that is given by, 
OT 


(264) 


T 4rcmyH 


Rewriting this in terms of the Eddington mass accretion 
rate (equation 241), 


M GM 


Mgaa ne (26) 


Ttrap = 
we see that photon trapping should become important 
roughly when the mass accretion rate first exceeds the 
value where the nominal luminosity is the Eddington 
limit. At higher m, the radiative efficiency will neces- 
sarily drop as photons produced close to the black hole 
are trapped and advected across the event horizon with- 
out any chance to escape to infinity. Although derived 
here in spherical symmetry, the same physics applies to 
geometrically thick disk flows (Abramowicz et al., 1988; 
Begelman and Meier, 1982). This regime can now be 
simulated using radiation hydrodynamics (Jiang et al., 
2019b). 

Stepping back from the details of the last few pages, 
Figure 19 summarizes what we have deduced about the 
qualitative structure of black hole accretion flows. The 
crucial parameter—at least as far as energetics go—is 
the accretion rate in units of the Eddington mass accre- 
tion rate. The thin disk solutions of Shakura and Sun- 
yaev (1973) and Novikov and Thorne (1973), although 


they are surely wrong in many details due to their ad 
hoc treatment of angular momentum transport, probably 
provide a semi-quantitative description of disk structure 
across about three orders of magnitude in accretion rate, 
10-3 < M/Mpaa Ê 1. At both extremities, the thin 
disk solutions fail and the expected outcome of accretion 
with non-zero angular momentum is a geometrically thick 
disk. The cause is inefficient radiative cooling, at low m 
because the low density keeps the energy locked away in 
the ions, which cannot radiate, and at high m because the 
high density traps photons within the inflowing gas. Out- 
flows are likely to accompany both geometrically thick 
disk regimes, which are surprisingly similar despite the 
different physics at low and high m. Net magnetic flux 
is an additional parameter that can qualitatively modify 
the properties of both thick and thin accretion disk flows. 
What happens if the accretion rate is pushed even fur- 
ther to the extremes? At extremely low m, the plasma is 
highly collisionless, and the long mean-free-path means 
that processes such as anisotropic conduction and vis- 
cosity are potentially important for the bulk flow (Chan- 
dra et al., 2015). At extremely high m, cooling due to 
neutrino emission opens up a second regime of at least 
moderately thin disk accretion (Chen and Beloborodoy, 
2007; Di Matteo et al., 2002; Popham et al., 1999). High 
enough accretion rates, exceeding M ~ 107? Mo s~! 
onto a stellar mass black hole, are realized during the 
collapse of the cores of massive stars to form a black hole 
and a massive disk (MacFadyen and Woosley, 1999). 


V. WAVES IN DISKS 


Consider the following situations, 


collisionless 
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FIG. 19 Illustration of the main regimes of black hole accretion as a function of the accretion rate, scaled to the accretion rate 
that would generate the Eddington limiting luminosity for radiatively efficient accretion. (There is a strong argument that this 
illustration should have a second dimension, with the extra parameter being the net magnetic flux that threads the disk. The 


regimes above are appropriate if the net flux is relatively weak.) 


e A planet, on a circular orbit, interacts gravitation- 
ally with the surrounding gas in the disk. 


e A binary black hole merges due to the emission 
of gravitational radiation, which leaves the system 
at the speed of light. From the point of view of a 
surrounding circumbinary disk, the merger leads to 
a near-instantaneous reduction in the central point 
mass from M to M(1—e), with € being of the order 
of a few percent. 


e A star-star flyby warps the outer part of a disk 
around one of them into a non-planar shape. 


Although physically quite different, these are all situa- 
tions where perturbations excite waves within an accre- 
tion disk. Disk waves assume a variety of forms, due to 
the different forces involved (pressure, rotation, gravity, 
magnetic fields) and because the perturbations of physi- 
cal interest differ from system to system. A book could be 
written on the mathematical description of these waves, 
and in fact at least one book has been written (Kato, 
2016). Here, we briefly describe some simple examples. 


A. Modified sound waves 


The simplest disk wave is an axisymmetric distur- 
bance in an unmagnetized, non-self-gravitating, two- 
dimensional fluid. We have already done the work needed 
to understand this situation. Dropping the term that re- 
sults from disk self-gravity, the dispersion relation (equa- 


tion 79) is, 


w? =k + Pk, (266) 
where w is the frequency of the wave, k is the wavenum- 
ber, cs is the sound speed, and «k is the epicyclic fre- 
quency, 


dQ 


k? = 40? + 2rO_., (267) 
dr 


For a Keplerian disk, x = Qg. This dispersion relation 
describes sound waves whose properties are modified by 
the Coriolis force. The modification takes the form of a 
cutoff—as we go to large wavelengths (small k) the wave 
frequency w + +k, with the opposing signs correspond- 
ing to inward and outward propagating waves. Provided 
that x? > 0 (as is assuredly the case for a near-Keplerian 
disk) the wave frequency w is always real, so there is no 
instability. 


B. The linear wave equation 


The properties of more general waves in disks can be 
derived in multiple ways. The derivation given in the 
review by Balbus (2003) is quite straightforward and di- 
rect, and we follow it here. The starting point is the fluid 
equations in the absence of magnetic fields or viscosity. 
Sticking to the inertial rather than shearing sheet frame, 


the continuity and momentum equations are, 


Op 
Æ +v- (pv) = 0, (268) 
wy Veen age (269) 


To keep matters (relatively) simple, we adopt a poly- 
tropic equation of state, 


P= Kp; (270) 
with K and y constants. Although we will not need it in 
order to look at wave properties, the hydrostatic equilib- 
rium structure implied by this equation of state is, 


1_ = page o 
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where po is the mid-plane density. Unlike in the isother- 
mal case, the disk has a free surface at finite z. 

To obtain an equation for linear waves within disks, 
the procedure is to first linearize the relevant equations, 
and then convert them from partial differential equa- 
tions to algebraic equations by specifying a form for the 
perturbations. We can illustrate the idea by working 
through the steps for the continuity equation in detail. 
Adopting cylindrical polar co-ordinates (r, ġ,z), we as- 
sume that there is some equilibrium background disk 
structure p(r,ġ, z), v(r,ġ,z) in which the fluid velocity 
is purely azimuthal (i.e. we ignore any inflow). On top 
of these fields we add perturbations, i.e. we take, 


(272) 
(273) 


p> pt op, 
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that are assumed to be small. For the density this is as 
simple as requiring that |dp|/p < 1. Substituting in the 
continuity equation, the terms involving the background 
state drop out and the one second order term in the small 
quantities is neglected, 
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(274) 
Allowing the density perturbation to have an arbitrary 
radial and vertical form, but requiring it to be periodic 
in azimuth, we write, 
dp = ôp(r, z) expli(md — wt)], (275) 


where w is the frequency of the wave and m is an integer. 
This implies that, 


(276) 


In the equilibrium state we have purely azimuthal flow 
v = (0,vg,0). Using the expression for the divergence 
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operator in cylindrical co-ordinates we find that, 


10 
V - (pv) = —=; (vg6p) = im “2 5p =imQdp, (277) 
r r 
where Q is the angular velocity (which is not necessarily 
Keplerian). Collecting the surviving terms together, the 
final result is, 


— iw@dp + V - (pôv) = 0, (278) 
where we have defined, 
@=w—moQ. (279) 


The quantity ©, which comes up often in this sort of 
analysis, is called the “Doppler-shifted wave frequency” . 

Turning to the equation of state (equation 270) we de- 
fine the enthalpy as, 
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where caq is the adiabatic sound speed. Differentiating, 
we obtain, 


ôH = ALA (281) 
p 
as the relation between perturbations in density and per- 
turbations in enthalpy for this equation of state. 
For the momentum equation, use of equation (67) gives 
the leading and first order terms that result from the 
convective operator in cylindrical co-ordinates as, 


[(v + dv) -V] (v + ôv) = 
Ug OU, vg Qugdug in avg 
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The linearized momentum equation then gives, 
06 
—iwdu, — 20dug = -2H (283) 
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where, 
dn? 
kK? = 407 + ; (286) 
dlnr 


defines the epicyclic frequency. Solving for dv, and dug 
between equation (283) and equation (284) we have, 
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Once the rotation profile of the background disk flow is 
specified (giving us Q and x”), these equations relate the 
velocity and enthalpy perturbations for waves of partic- 
ular frequency w and azimuthal wave number m. 

The final step in the derivation is to insert the ex- 
pressions for dv into the linearized continuity equation 
(equation 278). Using equation (281) the result simpli- 


fies to, 
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Here, we have defined, 


Der -0?. (291) 
This is the governing equation for linear hydrodynamic 
waves in a disk, ignoring magnetic fields and disk self- 
gravity!!. Even with these simplifications, it’s rather 
complicated! What is easy to see, however, is that there 
are two ways in which the denominators in the equa- 
tion can equal zero. These are resonant locations. When 
D = 0 we have the condition for Lindblad resonance, 
while © = 0 defines the location of corotation resonances. 
These resonances are central to the theory of planet-disk 
interactions, because for a disk perturbed by an exter- 
nal potential angular momentum transfer occurs only in 
their vicinity (Goldreich and Tremaine, 1979). 


C. Density and inertial waves 


Suppose for now that we are not in the vicinity of a res- 
onance. Then equation (290) describes freely propagat- 
ing waves within the disk. It can be solved numerically, 
but we can gain analytic insight with the aid of some 
additional approximations. Assume a WKB solution of 
the form, 


LS 
H = A(r, z) exp [5] , (292) 
€ 
where ¢€ is a small parameter that is introduced so that 
the phase S/e varies rapidly. The radial and vertical 
wavenumbers are, 


Os 
kr = (293) 
Os 


11 Note that I think there are typos in the signs of the last two 
terms in Balbus (2003). Check yourself before blindly using this 
equation. 
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FIG. 20 The dispersion relation for disk waves (equation 296). 
Plotted are contours of constant Doppler-shifted wave fre- 
quency @, in the (kz, kr) plane. The wavenumbers are scaled 
by Caa/K, where caa is the adiabatic sound speed and « is 
the epicyclic frequency. The blue ellipses describe high fre- 
quency density waves, with ©? > «?. The curves are shown 
for @/« between 1.1 and 2.0, in intervals of 0.1. The red hy- 
perbolae describe low frequency inertial waves, with ©? < kK’. 
The curves are shown for @/« between 0.1 and 0.9, again in 
intervals of 0.1. 


We further assume that k, > m/r, kz >> m/r, and that 
the disk is thin, so that caa < rQ. Substituting in equa- 
tion (290), and keeping only those derivative terms of 
order €~?, the leading order result is, 


Ap ; * Ap; a Apis 
LP iSe (5) P piS/e (F) — ee (295) 
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The dispersion relation is then, 
kica 4 KGa c], (296) 
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By taking the WKB limit the general wave equation re- 
duces to a pretty simple dispersion relation. 

Recall now some basic properties of conic sections. El- 
lipses are described by xz?/a? + y?/b? = 1, hyperbolae 
by x?/a? — y?/b? = 1. Depending upon the Doppler- 
shifted wave frequency, solutions to equation (296) in 
the (kr, kz) plane can be of either form. High-frequency 
waves, with ©? > «?, trace out ellipses in this plane, while 
low-frequency waves, with ©? < «Kĉ, trace out hyperbo- 
lae. These waves evidently have quite different geometric 
properties, and in fact represent two distinct classes of 
waves that can be present within disks. 

Writing equation (296) in the equivalent form, 


Dt — (keea +K’)? + 67k? c?, = 0, (297) 


where k? = k? + k?, the two types of disk waves can 
be defined by reference to which terms in the dispersion 
relation are dominant. Density waves are what we get if 
the third term in the dispersion relation can be neglected, 
either because k, = 0 (or is otherwise small) or because 
kcaa >> K. In this limit, 


D = k? + eak’, (298) 


and the physics is the same as §V.A. These are modified 
acoustic waves. Inertial waves, on the other hand, occur 
when the first term in the dispersion relation is negligible. 
Taking the limit where the sound speed cag —> œo, the 
dispersion relation takes the form, 


a2 k? Ke 
k2+k2 ` 
Inertial waves are low frequency disk phenomena, and 


unlike density waves involve fluid motions that are almost 
incompressible. 


(299) 


VI. WARPED AND ECCENTRIC DISKS 


It has been known for a very long time that not all 
accretion disks are circular, planar structures. The X- 
ray source Hercules X-1 exhibits a 35 dy periodicity that 
is interpreted as the precession period of a warped disk 
surrounding the accreting compact object (Katz, 1973). 
Maser emission from the nuclear region of the Seyfert 
galaxy NGC 4258 traces a geometrically thin and mildly 
warped disk down to sub-pc scales from the black hole 
(Miyoshi et al., 1995). More recently, several examples 
of warped (and in some cases “broken”) protoplanetary 
disks have been observed (e.g. Casassus et al., 2015; 
Kraus et al., 2020). Direct evidence for eccentric disks 
is sparser, but photometric variations known as “super- 
humps” in some cataclysmic variables imply excitation of 
the outer disk into an eccentric, precessing state (White- 
hurst, 1988). The formation of an eccentric disk also 
seems to be an inevitable consequence of tidal disruption 
events in galactic nuclei (e.g Bonnerot et al., 2016). 

The physics of warped and eccentric disks is subtle, 
and can be a viperous pit for the unwary. This Section 
gives a high level overview. 


A. Warped disks 


The state of a warped disk is jointly described by the 
surface density U(r) and the unit vector normal to the 
local disk plane Î(r)!?. We will assume that the disk is 


12 Some literature makes a distinction between warped disks, in 
which the ensemble of tilt vectors for the annuli in the disk lie 
in a single plane, and twisted disks, in which they do not. We 
won’t bother. 
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geometrically thin, that the orbits at all radii are circular 
(up to corrections due to the radial gas flow), and that 
I(r) is a smooth, continuous function. It is not immedi- 
ately obvious whether the evolution of the tilt vector, un- 
der the action of internal (and possibly external) torques, 
preserves continuity, but we defer for now discussion of 
the possibility that the disk breaks up into disjoint pieces 
(Nixon and King, 2012). 

Our immediate goal is to derive equations for the time 
evolution of the surface density and tilt of a warped disk. 
In addition to the regular planar disk physics there are 
new effects to consider, illustrated in Figure 21. If we 
consider two adjacent disk annuli, any mutual inclination 
results in vertical as well as radial shear. That vertical 
shear, in turn, leads to a periodic vertical displacement 
of the mid-planes of the two fluid annuli, which induces 
periodic radial pressure gradients. In a Keplerian poten- 
tial, the forcing is resonant with the epicyclic response of 
the disk, such that the misalignment launches a wave. 

What happens next depends upon the conditions 
within the disk. If the disk is sufficiently viscous, the 
nascent wave is damped locally, and we say that the 
warp evolution is in the viscous regime. The condition for 
viscous warp evolution is that (Papaloizou and Pringle, 
1983), 


a>, (300) 
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< h/r, radial communi- 
cation of the warp occurs via waves provided that the 
potential is nearly Keplerian, specifically that, 


In the opposite limit where a < 


Q? — 2 
eo < (301) 


h 
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The evolution equations for viscous and wave-like warped 
disks are different, and they predict qualitatively differ- 
ent behavior. Both regimes are physically relevant. Pro- 
toplanetary disks are thick and inviscid enough as to fall 
almost always into the wave-like regime. Disks around 
black holes and other compact objects, on the other hand, 
are expected to display a large variation in h/r with ra- 
dius (Figure 14). The outer regions (at a minimum) are 
likely to be described by viscous warp dynamics. 

Key references for the analytic theory of the fluid dy- 
namics of warped disks include Papaloizou and Pringle 
(1983), Papaloizou and Lin (1995), and Ogilvie (1999a). 
The review by Nixon and King (2016) provides a good 
overview, and is worth reading before wading into the 
technical details of the earlier papers. 


1. Viscous limit: heuristic derivation 


Conservation laws strongly constrain how warped disks 
behave. In the viscous limit, Papaloizou and Pringle 
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FIG. 21 Illustration (after Lodato and Pringle, 2007, their Figure 10) of the additional hydrodynamic effects present in a 
warped disk. Considering two neighboring but misaligned fluid annuli, there is obviously a vertical shear between them that is 
induced by the warp. As shown in the insets, the relative vertical displacements of the mid-planes of the annuli also result in 
time- and height-dependent radial pressure gradients. If the potential is close to Keplerian, the radial forcing is resonant with 
the disk’s epicyclic frequency. This physical effect means that a warped disk is a source of bending waves, which—depending 
upon the disk viscosity—may damp locally or propagate globally. 


(1983) and Pringle (1992) derived an evolution equation 
for geometrically thin warped disks by combining conser- 
vation laws for mass and angular momentum with an as- 
sumption that a distinct viscosity acts to damp misalign- 
ments. Although this approach does not quite recover the 
same equations as a full hydrodynamic analysis (Ogilvie, 
1999a; Papaloizou and Pringle, 1983), the differences are 
fairly minor, and the Pringle (1992) derivation is simple 
and physically instructive. We follow it closely here. 

We consider a geometrically thin disk with orbits that 
are circular but possibly misaligned with respect to each 
other. An annulus of the disk at radius r has width 
Ar, surface density X, radial velocity v,, and angular 
velocity Q. The tilt of the annulus in some inertial frame 
is described by a unit vector l that is normal to the disk 
(and parallel to the local angular momentum vector). 


Conservation of mass gives, 


A + a (rXv,) = 0. 


Ot | r Or (a02 


This equation is identical to the one for a planar disk. 

The angular momentum density (i.e. per unit area) is 
Er?QÎ, and the total angular momentum content of the 
annulus is then 2rrAr x Nr2Q1. The annulus’ angular 
momentum changes due to: 


e Mass flow into / out of the annulus. 
e Viscous torques acting within the disk. 


e External torques, present for example near Kerr 
black holes, around oblate central objects, and in 
binary systems. 


Mathematically (Papaloizou and Pringle, 1983), 


a (2nr’azÎAr) = 


2rr’ Qv, — Irr’ Q5 Înv, 


r—Ar/2 r+Ar/2 
+G (r + Ar/2) — G(r — Ar/2)+T. 


(303) 


G(r,t) is the viscous torque exerted on one annulus by 
its neighbor, T is the external torque. We can ignore T 
for now, as it can easily be added back in to the final 
evolution equation. 

Determining the form for G is the crux of the deriva- 
tion. In a warped disk there are two sources of shear, 
which we consider independently. First, there is an in- 
plane component, corresponding (in obvious notation) to 
the (r,@) stress, which acts in the direction of Î. This 
component has the same form as for a planar disk, 

Gi = See A. (304) 
dr 
Second, there is a component corresponding to the (r, z) 
stress, which has a cos ¢ azimuthal dependence, and van- 
ishes if 01/Or is zero. This component of the stress takes 
the form (Papaloizou and Pringle, 1983), 


G2 = 2nr(vo/2)Br05°r. (305) 
The factor of two comes from averaging the oscillatory 
vertical shear around the orbit. In this telling, vı and 
Və are separate viscosities that do not have a simple or 
obvious relationship. In particular, there is no reason 
to assume that vı = v2, and equality of these heuristi- 
cally defined “viscosities” is not what is expected for a 
fluid disk described by the Navier-Stokes equations with 


an isotropic viscosity. We will return to these possible 
pitfalls later. 

Using these forms for the viscous torque, and taking 
the limit as Ar > 0, the angular momentum conservation 
equation becomes, 


Ot r 


2 (=r?) + a (zuri) = 
1a (aze $i) pl? (zama) . (306) 
T 
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In addition to terms describing advection of angular mo- 
mentum (and thus tilt), the final term has a diffusive 
character with a diffusion co-efficient that is œ və. It ex- 
presses the tendency of viscosity to flatten out an initially 
warped disk. 

The combination of the vector equation (306) with the 
scalar equation (302) for the surface density completely 
determines the evolution of viscous warped disks, in the 
absence of external torques and given the specified as- 
sumptions. To combine them, v, has to be eliminated. 
We first dot Î into equation (306), and make use of the 
following identities which follow from the fact that Îisa 
unit vector, 


1-1=0, (307) 
pA i (308) 
i. ul o. (309) 
The result is, 
a (Sr?Q) + - (Evr*Q) = 
(a) + serra, a (310) 


We now multiply the continuity equation (302) by r°Q 
and subtract it to eliminate the time derivative from the 
left-hand-side. For the final term we note that, 


ð f; a 4 OF 
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so 1- 82Î/ðr? can be replaced with —|01/dr|2. This gets 
us to, 


2 


i 
di] _ 9 


= (311) 


(312) 


The left-hand-side looks as if it involves a spatial gradi- 
ent of v,, but on application of the chain rule it simplifies 
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considerably to give just v,©O(r?Q)/Or. The radial ve- 
locity is then, 


x12 
o (4, Er? S2) — inDr’Q ‘ 


ai 
Or ðr 
ð 
Dr Ər (r?Q) 


(313) 


Up = 


The novel piece here is the second term involving vg and 
the radial gradient of the tilt vector. This term reflects 
the fact that viscous flattening of local variations in the 
disk tilt involves energy dissipation, which must be ac- 
companied by radial inflow of gas through the disk. Using 
this expression for vp, we can obtain separate equations 
for the time evolution of the surface density and tilt vec- 
tor. These can be found in Pringle (1992). At this point, 
physically, we are done. Noting that 1 is a unit vector, 
however, we still have four equations for only three in- 
dependent quantities. Some additional simplification is 
therefore possible. Specializing to a Keplerian potential, 
straightforward manipulations yield an evolution equa- 
tion for, 


L = (GMr)? si, (314) 
which takes the form, 
Bros a 
-_ emit (315) 


If the disk is planar, this equation reduces to the usual 
evolution equation for the surface density of a thin disk. 


2. Viscous evolution equation 


Imagine a minimal planetary system made up of two 
planets on stable, circular, but mutually inclined orbits. 
Gravitational torques will result in precession of the or- 
bital planes, but there is no dissipation in the system 
and both the energy and the angular momentum stay 
fixed. Replace the planets with interacting fluid rings, 
and it becomes obvious that elementary considerations of 
mass and angular momentum conservation cannot fully 
determine how a warped disk evolves, because such con- 
siderations cannot capture possible precessional effects. 
A comprehensive description must instead be based on 
analysis of the Navier-Stokes equations!®. A program to 


13 Assuming, as is customary in almost all analytic work, that we 
can’t provide a better simple description of angular momentum 
transport in a fluid disk. 


do so was started by Papaloizou and Pringle (1983) and 
essentially completed by Ogilvie (1999a). It leads to both 
a modified version of equation (315) and, more impor- 
tantly, to an understanding of the relationship between 
vı and 1. 

Ogilvie (1999a) is stiff mathematical medicine, which 
the reader seeking full details will have no choice but 
to imbibe. A central result is the modified version of 
equation (315). Keeping the notation the same, it reads 
(Nixon and King, 2016), 
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dt rdr| = ðr 
o Die a L 
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As presaged above, the modification takes the form of an 
additional term which describes precession driven by the 
presence of a warp. Its magnitude depends on a function 
v3, which has the dimensions of a viscosity. In some cases 
this extra term is negligible, but there is no reason (other 
than a rather minor simplification) not to include it in 
calculations of the evolution of warped, viscous disks. 

The relation between the dissipative vertical and hor- 
izontal viscosities, v2 and vı, can be determined from 
the hydrodynamic theory. For small amplitude (lin- 
ear) warps the result is (Ogilvie, 1999a; Papaloizou and 
Pringle, 1983), 


v2 2(1+7a7) 1 
vı at(4+a2) 2a?’ 


(317) 


where the final expression is valid for the usual case where 
a & 1. It is worth emphasizing that this normally large 
ratio between vz and vı occurs even for an intrinsically 
isotropic fluid viscosity, and has nothing to do with any 
anisotropy that might arise, for example, if angular mo- 
mentum transport originates with the MRI. 

The local strength of a warp in an accretion disk de- 
pends upon the dimensionless quantity, 


al 


Pak (318) 


v= 


The analysis of Ogilvie (1999a) yields equations that can 
be solved to obtain the viscosities (11, v2 and v3, or Qı 
through Q3 in Ogilvie’s notation) as f(w). They have a 
non-trivial form, with both vı and v2 typically decreasing 
as the warp amplitude becomes increasingly non-linear. 
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3. Bardeen-Petterson effect 


Returning to equation (315), we now consider the effect 
of external torques T on the evolution of warped disks. 
The most famous example is the Lense-Thirring effect 
(Lense and Thirring, 1918), which leads to the precession 
of orbits around Kerr black holes that are inclined with 
respect to the equatorial plane. The precession frequency 
can be written, approximately, as, 


(319) 


where wyr = 2GJ/c?, and J = aGM?/c is the angular 
momentum of a black hole with dimensionless spin pa- 
rameter a. Bardeen and Petterson (1975) observed that 
differential precession due to the Lense-Thirring effect, 
acting on an initially planar but tilted disk, would rapidly 
warp the disk close to the black hole and drive up |W}. 
Viscosity would then damp the disk into the equatorial 
plane at small orbital radii, while the disk further out 
remained misaligned. Provided that the disk is in the 
viscous limit, the same sort of effect can occur when disk 
precession is driven at small radii by the oblateness of a 
central body, or at large radii by torque from a binary 
companion. 

Starting from the equations of Papaloizou and Pringle 
(1983), Kumar and Pringle (1985) computed the shape 
of a warped disk subject to the Bardeen-Petterson effect. 
A simple analytic solution for the shape can be derived 
in a limit where the warp is small, the viscosity is con- 
stant, and the structure of the disk close to the inner 
boundary is not important (Scheuer and Feiler, 1996)"4. 
The starting point is equation (315). We drop the second 
order term |01/dr|?, set the time derivative to zero, and 
add in a term that represents forced precession due to 
the Lense-Thirring effect. In steady-state, a warped disk 
around a spinning black hole obeys, in the viscous limit, 
the equation, 
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= rdr 
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(320) 


For now, we pick co-ordinates such that the spin of the 
black hole is aligned with the z-axis, so that the vector 
wr in the above equation is, 


wur = wrr (0,0,1), (321) 
2GJ GM\? 
wr =a = 2ac (S) (322) 


14 Peter Scheuer taught one of my undergraduate mathematics 
courses, and this argument is a nice example of a model mathe- 
matical physics problem. 


Task one is to determine the surface density profile of the 
disk. For a warped disk this is expressed via the angular 
momentum density, 


L = Li, 
L= |L| = (GMr)"? Ð. 


(323) 
(324) 


In terms of these variables, equation (322) becomes, 


1d 3r d 3 
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(325) 


To convert this to a scalar equation for L we use the same 
trick as before. We dot through with 1 and make use of 
the fact that 1- dl/dr = 0. The result is, 


3 
3r— (vı L) Sn ; (326) 
Assuming that vı is independent of radius, the solution 
is, 


L= cor"? — 2c, (327) 
where co = (GM)!?E» and Doo is the surface density 
at large disk radii. cı is a constant that depends on 
the inner disk boundary condition. The result, that for a 
constant viscosity X tends to a constant as r — oo, is the 
same as for a planar disk (equation 119). Small warps do 
not affect the steady state run of the disk surface density. 

We now turn to the crux of the problem: determining 
the shape of the warp. Using the solution for L (equa- 
tion 327) in equation (325) gives, 


d ‘4 1/2 aĵ 
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(328) 
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where we have again assumed that vı is not a function 
of radius. Additionally, assume that the disk region of 
interest is sufficiently far away from the black hole that 
the inner disk boundary condition is immaterial. This 
allows us to drop cy, leaving, 


d [1 gdl] 
dr jem 1|- 
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(329) 


Writing 1 = (ls, ly, l), wur x Î = wur(—ly,lzx,0). Be- 
cause Îis a unit vector, only two components are indepen- 
dent, and it is convenient to cast the equation in terms 
of a complex variable, 


W = l, + ily. (330) 
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FIG. 22 The shape of a warped disk around a Kerr black 
hole, computed in the viscous limit (with radially constant 
viscosities) using the Scheuer and Feiler (1996) solution. The 
blue, red and black curves show the components of the tilt, 
lz and ly, and the total tilt 1 = ,//2 + l2. For this example, 
(wir /V2rin) = 300. The Scheuer and Feiler (1996) solution 
applies to small warps, and hence the values on the y-axis are 
arbitrary. An accretion disk subject to the Bardeen-Petterson 
effect is predicted to be warped across a broad radial range, 
and to be modestly twisted. 


We obtain, for constant 12, 
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This can simplified further with the substitution X = 
2r—"/2 to give, 
d? W 2iwyr W 
dx? — V2 ` 


(332) 


The solution is (Scheuer and Feiler, 1996), 


W = kexp f (i—1) (=) , (333) 


ar 


with k a constant. As r => co we have W — k, represent- 
ing a flat inclined disk, while W — 0 at small radius. The 
solution is well-behaved even though, as we have noted, 
we have dropped terms that are non-negligible close to 
the black hole. The solution is plotted as Figure 22. 
The shape of a disk subject to the Bardeen-Petterson 
effect is not something that is easily observed, so a fair 
bit of the interest in the problem lies in what the warped 
disk does to the black hole. Given enough time, we expect 
the angular momentum vectors of the black hole and disk 
to change, until we reach an end-state with a flat disk 
in the equatorial plane defined by the (final) black hole 
spin vector. The disk may be either aligned or counter- 
aligned with the black hole. General considerations of 
angular momentum conservation provide constraints on 
which of these outcomes occurs (King et al., 2005), but 


to get at the time scale for realignment a disk model is 
indispensable. We sketch out how it’s done in the case 
of the Scheuer and Feiler (1996) solution. For the disk, 
with angular momentum Ją, the rate of change of angular 
momentum due to the Lense-Thirring effect is, 


Ta n faxt 
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= 2TiC2WLT Í Wr-9/*dr. (336) 
Substituting for W we can do the integral, with the result 
being, 
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(337) 
We now consider the black hole’s spin. Its angular mo- 
mentum changes at an equal and opposite rate to that of 
the disk, 

dJpu 


Ia 
dt — dt’ 


(338) 


Switching to a disk-centered coordinate system, the disk 
as r — œ lies in the x-y plane, rather than that plane 
coinciding with the equator of the black hole as before. 
Then k = —(jz + ijy), where j, and jy are components 
of the unit vector describing the black hole spin. Equa- 
tion (337) then gives the time evolution of the black hole 
spin as (Scheuer and Feiler, 1996), 


d(e + ijy) _ 
dt 
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The imaginary part of this equation describes precession 
of the black hole, while alignment of the black hole spin 
with the angular momentum of the disk at large radii is 
described by the real part. Both precession and align- 
ment occur on the same time scale, 


1 facm\'? 
T doo ( Gi ) 
The dependence on disk parameters comes in via v2 and 
Xo (which, for a given steady accretion rate M, would be 
inversely proportional to 1). 

This illustrative version of the analysis can be extended 
to account for radially varying viscosities (Martin et al., 
2007; Natarajan and Armitage, 1999). Gerosa et al. 
(2020), who also study the effect of precession due to a 
binary companion, provide a detailed analysis. The nu- 
merical value of the alignment time scale (equation 340) 
depends on the details of the disk model, but generically 
it works out to be short—compared in particular to the 
Salpeter time scale on which accretion would change the 
magnitude as well as the direction of J. The ease with 


tp = talign = (340) 
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which a misaligned disk can change the spin axis of a 
black hole occurs because of the larger specific angular 
momentum of disk gas, and because in the viscous regime 
we expect v2 > vı. Natarajan and Pringle (1998), for 
example, estimate that talign may be of the order of one 
percent of the Salpeter time. 


4. Wave-like evolution equation 


Warps in protostellar disks, and in the inner regions of 
AGN disks, are expected (according to equation 300) to 
evolve in the wave-like regime, and the prior viscous anal- 
ysis no longer holds. Viewed as m = 1 bending waves, 
the linear evolution of warps is a subset of the general 
problem of disk waves (§V). It has been treated by Pa- 
paloizou and Lin (1995), Demianski and Ivanov (1997), 
and Lubow and Ogilvie (2000). These papers use differ- 
ent notation and address distinct scientific questions, but 
fundamentally solve the same linear problem. As before, 
Nixon and King (2016) is recommended for an accessible 
introduction. 

The derivation of the equations for the linear evolution 
of warp waves in an accretion disk is given, as compactly 
as possible, as an Appendix in Lubow and Ogilvie (2000). 
We won’t repeat it here. If G(r,t) is the internal torque 
in the disk, the evolution of the unit tilt vector 1 in the 
absence of external torques is given by the coupled equa- 
tions, 


al 10G 
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The right-hand-side of the second of these equations has 
some dependence on the disk’s vertical structure (Lubow 
and Ogilvie, 2000), with the form given here being true 
if it is isothermal. For a strictly inviscid disk (a = 0) we 
can combine these equations into a single wave equation 
for the disk tilt, 
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With a rescaling of the radial co-ordinate (e.g Ogilvie, 
2006) this equation can be cast into the form of the clas- 
sical wave equation, and one finds that linear bending 
waves propagate at speed (Papaloizou and Lin, 1995), 
Cs 
warp — >° 344 
Y P 2 ( ) 


These waves are non-dispersive. 


5. Additional dynamical considerations 


Several other dynamical or hydrodynamical effects can 
come into play for warped disks. We list a few, with brief 


comments, here. 


The dichotomy of having, separately, a non-linear the- 
ory for viscous disks with a > h/r, and a linear theory 
(only) for wave-like warp evolution with a < h/r, is ob- 
viously unsatisfactory. How does a disk with a = h/r, or 
one where different radial parts fall into different regimes, 
evolve? Martin et al. (2019) proposed an evolution equa- 
tion that unifies the two regimes, valid in the limit where 
the warp is small. Their work was extended by Dulle- 
mond et al. (2021), who used a formalism developed by 
Ogilvie and Latter (2013) to derive a closely related equa- 
tion describing warped disk evolution. 


We introduced this section by noting that the key phys- 
ical reason why warped disks behave differently from pla- 
nar ones is because the warp drives oscillatory, and of- 
ten resonant, radial gas flows (Figure 21). That these 
secondary flows can be unstable was suggested by Pa- 
paloizou and Terquem (1995) and studied in detail by 
Gammie et al. (2000). Deng et al. (2021) have sim- 
ulated the resulting parametric instability in a global 
context, finding that the resulting turbulence strongly 
damps warps. 


A warped disk is “held together” by internal torques, 
which may not always be strong enough to prevent the 
disk breaking up into disjoint annuli. Disk breaking 
(Nixon and King, 2012) occurs whenever a disk is subject 
to particularly strong differential precession, for example 
due to the Lense-Thirring effect (Liska et al., 2021; Nixon 
et al., 2012) or due to torques from a binary companion 
(Nixon et al., 2013). It may also be possible for a warped 
disk, with an initially smooth radial profile of 1, to evolve 
slowly (e.g. on a viscous time scale) up to the point where 
regularity of the solution is lost, and a break develops 
(Dogan et al., 2018). The existence of this latter channel 
for breaking depends upon how the internal torques, vı 
through v3 in equation (316), vary with the local warp 
strength |y]. 


Warped disks are commonly found in binaries. The 
gravitational potential of a binary gives rise to rich purely 
gravitational dynamics, most famously the Kozai-Lidov 
effect (Kozai, 1962; Lidov, 1962), which leads to large- 
amplitude oscillations in eccentricity and inclination for 
test particles that are sufficiently inclined to the binary 
plane (for a review see Naoz, 2016). Martin et al. (2014) 
identified analogous dynamics in simulations of initially 
misaligned fluid disks in binary systems. Linear analysis 
(Lubow and Ogilvie, 2017; Zanazzi and Lai, 2017) shows 
that the fluid version of the Kozai-Lidov effect depends 
upon h/r, and that unlike the free particle version can 
be present even for low values of the initial misalignment 
angle. 
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B. Eccentric disks 


Returning to planar disks, we drop the prior assump- 
tion that orbits in the disk are circular (up to a small 
correction due to radial inflow) and instead allow them 
to be eccentric. This leads to new complications: 


e The time scale for establishing vertical hydrostatic 
equilibrium in an accretion disk, thyaro ~ 27! 
(equation 164), is of the same order as the orbital 
time scale. In an eccentric disk, the vertical ac- 
celeration due to the gravity of the central object 
changes on this time scale between apocenter and 
pericenter, dramatically so if e is large. Eccen- 
tric disks are therefore never in true vertical hydro- 
static equilibrium (Ogilvie and Barker, 2014), and 
the vertical structure couples strongly to the dy- 
namics in situations—such as disks resulting from 
Tidal Disruption Events—where the eccentricity is 
large (Lynch and Ogilvie, 2021la; Ryu et al., 2021; 
Zanazzi and Ogilvie, 2020). 


e Depending upon the nature of angular momentum 
transport within the disk, an initially circular disk 
may exhibit viscous overstability (Lyubarskij et al., 
1994; Ogilvie, 2001). If this is the case, the intu- 
itive expectation that an eccentric disk ought to 
eventually circularize would be violated. Exten- 
sion of the viscous a-model (with a Navier-Stokes 
viscosity) to eccentric disks is not recommended, 
and instead we should think about how magnetic 
fields and the magnetorotational instability inter- 
act with disk eccentricity (Chan et al., 2018; Lynch 
and Ogilvie, 2021b). A direct comparison of two- 
dimensional a disks with three-dimensional MHD 
simulations, in binary systems with mass ratios in 
the range where eccentricity is resonantly excited, 
shows that MHD effects work to damp eccentricity 
(Oyang et al., 2021). 


e Eccentric disks (in common with warped disks) ex- 
hibit a hydrodynamic parametric instability (Pa- 
paloizou, 2005; Wienkers and Ogilvie, 2018). This 
can generate turbulence, for example in circumbi- 
nary disks whose eccentricity is maintained by ex- 
ternal forcing (Pierens et al., 2020). 


Figure 23 illustrates example streamlines for fluid in an 
eccentric disk. Both the eccentricity e and the longitude 
of pericenter w may vary as functions of some radial co- 
ordinate used to label streamlines (the semilatus rectum, 
à = a(l — e?), where a is the semi-major axis of the 
orbit, is a good choice). These represent two additional 
functions whose evolution must be specified, unlike in the 
case of a thin circular disk where we have only the surface 
density. 

Although it does not address all of the complexities 
of eccentric disks, the linear theory of disk eccentricity 
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FIG. 23 Example streamlines for two eccentric disks. In the left-hand example, the eccentricity e varies radially, but the 
argument of pericenter w remains constant. In the more general case, shown on the right-hand side, both e and w are functions 
of distance from the central object. It is often convenient to use the semilatus rectum À = a(1 — e°) as the radial co-ordinate 


when describing eccentric disks. 


is already useful. It is the foundation of models for how 
planets and binary systems affect their surrounding disks. 
We sketch here the main steps of the derivation by Good- 
child and Ogilvie (2006). An incomplete list of related 
perturbative analyses include those by Kato (1983), Lee 
and Goodman (1999), Tremaine (2001) and Papaloizou 
(2002). 

Following Goodchild and Ogilvie (2006) we consider a 
two-dimensional inviscid disk model in cylindrical polar 
coordinates, with density p(r,@), velocity v = (vr, vo), 
pressure p, and gravitational potential ®. The poten- 
tial is assumed to be that of a central object, and to be 
axisymmetric!>. The fluid equations are, 
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15 Goodchild and Ogilvie (2006) include an additional term repre- 
senting a tidal potential, which we ignore here in the interests of 
thinking about the eccentricity evolution of isolated disks only. 


The pressure is related to the density with an adiabatic 
exponent y. 

The equilibrium state is a circular disk, with v, = 
0/0¢ = ð/ðt = 0. This gives the usual radial force 
balance, 


t= yp, (349) 


We now assume, here and subsequently, that the disk is 
thin such that cs < vg. The pressure gradient term is 
then negligible and, 


_ de 


= — 
$ dr’ 


(350) 
with the angular velocity Q = vgo/r being that in the 
unperturbed disk state. 

Our goal is to reduce equations (345) through (348) 
into a PDE for the time evolution of linear eccentricity 
perturbations to the base state. To that end, we pursue a 
variant of the sort of stability analysis we’ve done before, 
starting by linearizing and perturbing the fluid equations. 
We then Fourier analyze in the azimuthal co-ordinate 
only, restricting attention to the m = 1 mode that we 
will identify with eccentricity. Proceeding informally we 
make the substitutions, 


Ur > v, exp[—4d], ( 

vg > rQ + v, exp[—id], (352 
p> p+ p exp|~iġ], ( 
p — p + p' exp[-i] ( 


where the perturbed quantities denoted with primes are 
understood to be small. Ditching second order terms in 
the perturbations, 
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These manipulations have eliminated the azimuthal 
derivatives, and give equations for the perturbed vari- 
ables that are just functions of radius and time. 

Next, we look for the lowest order solution to the per- 
turbation equations. It’s fine to assume that the time 
scale for eccentricity evolution is long compared to the 
orbital time scale, so that |Ov/./Ot| < |iQu}|. The pres- 
sure terms are also negligible at lowest order for a thin 
disk, so the first two of the equations above both imply, 


2v = —iv,. (359) 
Accordingly, we can write, 


vl. = irQ E(r,t), (360) 


v, = Fre, t). (361) 
The perturbations to the radial and azimuthal velocity 
are out of phase, and differ in magnitude by a factor 
of two. For eccentricities e « 1 these characteristics 
describe a two-body eccentric orbit, and we can therefore 
identify E(r,t) with a complex eccentricity, 


E(r,t) = eexplia]. (362) 


Both the eccentricity e and the longitude of pericenter w 
are functions of radius and of time. 

At this point physically we are done. The perturbation 
equations can be combined to give a single equation for 
the time evolution of the complex eccentricity, 
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(363) 


The first term on the right-hand-side describes preces- 
sion in the presence of a radial pressure gradient, while 
the second describes radial diffusion of the complex ec- 
centricity. We can look for normal mode solutions to this 
equation by writing, 


E(r,t) = E(r) expliv, (364) 
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where w is an eigenvalue determined as the solution of, 
E dp 1d 3 dE 
+ ar (PG) 


2rQwE = ade 
In this highly simplified and inviscid model, the eccen- 
tricity evolves as a superposition of precessing normal 
modes, whose structure and amplitude are set by the ini- 
tial and boundary conditions. 

The reader who wants to delve further into the linear 
theory of disk eccentricity can find it in Goodchild and 
Ogilvie (2006), which contains a lot of interesting stuff. 
It turns out that including the effects of viscosity leads 
to only a small modification of equation (363), 
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Here ap is a Shakura-Sunyaev style bulk viscosity. It acts 
straightforwardly in the context of this theory to damp 
disk eccentricity. Stepping back, however, one might ask 
whether a Navier-Stokes viscosity (with shear and bulk 
terms) really applies to a turbulent accretion disk, and 
if not what sort of effective bulk viscosity is generated 
by the MRI or some other physical angular momentum 
transport process. These are assuredly not straightfor- 
ward questions. 


(365) 
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VII. CLASSICAL DISK INSTABILITIES 


Thus far, the instabilities we have discussed have been 
linear instabilities of the fluid equations (mostly in the in- 
viscid Euler limit), with the addition of magnetohydrody- 
namics or self-gravity. They are “disk” instabilities inas- 
much as the choice of unperturbed state is appropriate 
for Keplerian disk flow. Interest in their non-linear evo- 
lution centers on their role in angular momentum trans- 
port, turbulence, and diffusion. 

These fluid instabilities are not the only class of in- 
stability of interest. There is another. The equations 
describing disk evolution (for example the diffusive equa- 
tion for the surface density, equation 114), can also ex- 
hibit instabilities. These instabilities are less fundamen- 
tal than, say, the MRI, because the equations themselves 
have baked in various assumptions that are only approx- 
imate. They are nonetheless of great interest, first be- 
cause they underly what was historically a great success 
of accretion disk theory—the identification of thermal 
disk instability with observed dwarf nova outbursts—and 
second because they allow a simplified analysis of com- 
plex physical problems such disk warping. 


A. Viscous instability 


The condition for viscous stability is determined by 
first considering a steady-state solution (1) to the one- 
dimensional disk evolution equation (114). Following 


Pringle (1981) we make the substitution 4 = v% and 
consider perturbations u — u + ôu. We assume that 
v = V(X). Substituting in the evolution equation (114), 
the perturbation dp evolves as, 


Sy) = HBO eB (5). 
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The perturbation ĝu grows if the diffusion coefficient, 
proportional to ðu /Ə£, is negative. This defines the con- 
dition for viscous instability, which takes the simple form, 


a (vu) < 0. 


5 (368) 


Toy disk models where v is a constant, or a fixed func- 
tion of radius, are evidently stable, as would be expected 
given the diffusive nature of the governing equation. In 
principle, however, there is no reason why some physical 
mechanism for angular transport could not yield an effec- 
tive viscosity law that implied viscous instability. Such 
a disk would have an intrinsic tendency to break up into 
rings. How viscous instability would saturate is not en- 
tirely clear, though the sharpness of ring-like structures 
in disks is known to be limited by the onset of the Rossby 
Wave Instability (RWI; Lovelace et al., 1999)!6. The 
RWI, in turn, leads to the formation of non-axisymmetric 
structure in the form of vortices (Li et al., 2001; Richard 
et al., 2013). 


B. Thermal instability 


For a geometrically thin disk in thermal equilibrium, 


Q4 = Q-, 


where Q+ is the heating rate per unit surface area of the 
disk and Q- is the corresponding cooling rate. Consider 
a perturbation to the mid-plane temperature Tp, on some 
scale A >> h, so that cooling of the perturbation occurs 
predominantly vertically rather than radially. Thermal 
instability occurs if, 


dlog Q+ 
dlog Te 


(369) 


dlog Q- 
dlog Ts ` 


(370) 


When the condition for thermal instability is satisfied, 
either upward or downward perturbations to T, lead to 
runaway. 

The instability condition is readily evaluated for var- 
ious flavors of a-model disks. Starting with the right- 
hand-side, which represents cooling, Q- = Volos Using 


16 RWI physics was discussed earlier in a galactic context, under 
the monikor of the “negative mass instability” (Lovelace and 
Hohlfeld, 1978). 
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the one-zone relation between the central and effective 
temperature (equation 181), 
Tê T4 
Q- = 20Th x aa re (371) 
with « the (mid-plane) opacity. For a constant opacity, 
for example that due to electron scattering, we find that 
dlog Q_/dlog T, = 4. 
What about the left-hand-side? Using equation (120) 
together with the a-prescription, 


9 
Q4 = aÊ Dx. 


3 (372) 


For a gas pressure dominated fluid c? « T., and we 
find that dlogQ,/dlogT, = 1. The simplest limit 
of a gas pressure dominated disk, with a temperature- 
independent opacity, is thus thermally stable. 

Everything changes if radiation pressure is dominant. 
For a radiation pressure dominated gas, the sound speed 
2 x T!/p, and 


Q4 x Tth. (373) 


Noting that the scale height h scales with the cen- 
tral pressure, and hence with T+, we find that 
dlogQ+/dlogTe = 8. The inner, radiation pressure 
dominated / electron scattering opacity zone of the 
Shakura-Sunyaev model, is analytically expected to be 
thermally unstable! The presence of such an instability 
would be expected to drive large-amplitude short-time 
scale variability of accreting black hole systems (Light- 
man and Eardley, 1974; Shakura and Sunyaev, 1976). Al- 
though accreting black holes are indeed variable sources, 
the character of the observed variability does not match 
simple expectations for a thermal instability origin. In 
particular, stellar-mass black holes, whose accretion rates 
cycle between values where the inner disk would be either 
gas or radiation pressure dominated, do not show obvious 
changes in their variability properties at the point where 
thermal instability would be predicted to set in. 

The existence of thermal instability in radiation pres- 
sure dominated a disks depends upon multiple assump- 
tions hardwired into the a model, many of which are 
known to be false, at least when examined closely enough. 
Agol et al. (2001), and many subsequent authors, stud- 
ied the thermal stability of radiation pressure dominated 
disks using either viscous or MHD simulations. Current 
results suggest that weakly magnetized disks would ex- 
hibit thermal instability (Jiang et al., 2013; Mishra et al., 
2016). Disks threaded by enough net flux as to be at least 
marginally magnetic pressure dominated, however, are 
stable (Begelman and Pringle, 2007; Jiang et al., 2019a; 
Oda et al., 2009; Sadowski, 2016), and this is a possible 
resolution of the apparent discrepancy between observa- 
tions and theory highlighted above. 

The radiation pressure version of thermal instability, 
relevant for very hot disks, is driven by the behavior of 


the heating term in the thermal equilibrium equation. 
One can also get instability, at much lower tempera- 
tures, from the cooling term. From equation (371) it is 
clear that one can get a small value of dlog Q- /dlog Te, 
and thermal instability, if the opacity « is a sufficiently 
strongly increasing function of temperature. Physically, 
this occurs when hydrogen is partially ionized at Te ~ 
104 K, and the opacity is dominated by H~. Analytic 
fits to the opacity in this regime yield x « T° (e.g. Bell 
and Lin, 1994), which is easily strong enough to source 
thermal instability under gas pressure dominated disk 
conditions. 


1. The S-curve and limit cycles 


Thermal instability, associated with the ionization of 
hydrogen, is well-accepted as the root cause of outbursts 
in dwarf nova systems. Dwarf novae are a subclass of 
cataclysmic variables (Warner, 1995) in which a weakly 
or non-magnetized white dwarf accretes from a low-mass 
main sequence star in a Roche-lobe filling mass transfer 
binary system. The prototypical dwarf nova outburst, 
seen in systems such as U Geminorum, is observed as 
roughly 5 magnitude increases in optical brightness, that 
last about a week and recur every few months. 

The identification of dwarf nova outbursts with ac- 
cretion disk thermal instabilities was made in the 1970s 
(Hoshi, 1979; Osaki, 1974), and developed into a working 
quantitative model by many authors in the early 1980s 
(Cannizzo et al., 1982; Faulkner et al., 1983; Meyer and 
Meyer-Hofmeister, 1981; Mineshige and Osaki, 1983b; 
Smak, 1984). Excellent reviews include those by Lasota 
(2001) and Hameury (2020). 

Dwarf nova outbursts involve both local and global 
disk processes. At a local level, a thermally unstable 
annulus of the disk has an unstable thermal equilibrium 
for some given surface density © and angular velocity Ox. 
Such an annulus would evolve on the thermal time scale 
~ 1/(aQx) to one of two stable states: a cold state where 
the central temperature is such that hydrogen is substan- 
tially neutral, or a hot state where it is substantially ion- 
ized. We can visualize this state of affairs by plotting the 
thermal equilibria for disk annuli in the (£, Teg) plane. 
Figure 24, which is a schematic based on calculations 
by Bollimpalli et al. (2018), shows what this looks like. 
The thermal equilibria have a characteristic “S-curve” 
appearance, such that two stable equilibria (and an in- 
termediate unstable one) exist across a range of surface 
densities. The computed S-curves are a function of white 
dwarf mass, disk radius, and assumed a value(s). 

At this point it’s worth highlighting an important sub- 
tlety of the dwarf nova disk instability model. The ex- 
istence of a classical thermal instability, here due to the 
behavior of the opacity near the hydrogen ionization tem- 
perature, is not a sufficient condition for producing a 
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clear S-curve. Almost equally important are changes to 
the disk vertical structure that occur due to vertical en- 
ergy transport by convection. Even then, the computed 
S-curves are rather lackluster. Prominent S-curves, of the 
form shown in Figure 24, occur if we make the further as- 
sumption that a on the hot branch is significantly larger 
than on the cold branch. Assumed ratios Qhot/&cola © 10 
do the trick. 

Figure 24 shows qualitatively the limit cycle that devel- 
ops in dwarf novae disks as a consequence of the S-curve. 
During quiescence, gas from the Roche lobe filling sec- 
ondary streams onto the outer accretion disk at a rate 
that exceeds what the disk can transport while on the 
cold branch of the S-curve. Mass accumulates in the 
disk, until at some radius the surface density exceeds the 
maximum surface density, “max, permitted on the cold 
branch. Usually, but not always, the critical point is close 
to the white dwarf. Once © > “max the only allowable 
local thermal equilibrium solution for the disk lies on the 
hot branch. The annulus that has been triggered rapidly 
heats up on the thermal time scale ~ 1/(aQx). Radial 
diffusion of mass and heat can then act to trigger the 
cold — hot transition in neighboring annuli, and a heat- 
ing front sweeps through the disk until every annulus is 
on the hot branch. Once fully in outburst, accretion onto 
the white dwarf exceeds the rate of mass supply from the 
secondary, and the surface density drops, triggering a re- 
turn to quiescence. 

Quantitative models for disk instability rely on sep- 
arate calculations for the vertical structure and radial 
evolution (“1+1D”). In the radial direction, it is neces- 
sary to treat departures from thermal equilibrium, and 
therefore the usual diffusive equation for the surface den- 
sity (equation 114) is supplemented with an equation for 
the evolution of the central temperature Te, 


oTe Q+-Q-+I RT.18(rvr) OT. 
Oot Cp g 


[Cp T Or " Or (ara) 
Here R is the gas constant, u is the mean molecular 
weight, and c, is the specific heat capacity at constant 
pressure. J is the radial flux of energy due to basically 
diffusive processes, which could be radiative diffusion and 
/ or turbulent radial heat transport. Different forms 
for J have been used in the literature (Cannizzo, 1993; 
Hameury, 2020). 

The agreement between disk models and observations 
of dwarf novae (e.g. Cannizzo, 1993) is impressive”, but 
leaves open questions such as why the efficiency of an- 
gular momentum transport on the hot branch should 
be much higher than on the cold branch. Local nu- 
merical simulations suggest that convection enhances the 


17 After pages and pages of beating up on the œ model, you might 
have been wondering why anyone takes it seriously. This is a 
large part of the reason why. 
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FIG. 24 Schematic illustration of how the existence of a local S-curve leads to global outbursts in dwarf novae accretion disks. 
The S-curves shown for different disk radii in the X-T plane are simplified from calculations by Bollimpalli et al. (2018), 
for a 1.35 Mo white dwarf. The depicted global evolution, in which outbursts start from the inside-out and decay from the 
outside-in, is typical for dwarf nova models but not a guaranteed property. 


strength of MRI transport near the tip of the hot branch, 
while Ohmic diffusion can damp MHD turbulence on the 
cold branch (Coleman et al., 2016, 2018; Hirose et al., 
2014; Scepi et al., 2018). As in other accreting systems, 
net magnetic fields and disk winds may contribute to the 
observed behavior (Scepi et al., 2019). 


C. Radiation warping instability 


A warped disk that intercepts and re-radiates radia- 
tion from a central source experiences a radiative torque 
that modifies the warp. Pringle (1996) showed that the 
simplest model of this physical situation—a razor-thin 
viscous disk illuminated by a small central point source 
of radiation—is subject to a linear instability that can 
warp an initially almost planar disk. The derivation in 
Pringle (1996) is clear, and we won’t repeat it here. Fur- 
ther linear analysis was given by Maloney et al. (1996). 

The radiation warping instability sets in beyond a crit- 
ical radius, which depends upon the luminosity and vis- 
cosity of the disk. For a steady, self-luminous disk with 
radiative efficiency ņ, the critical radius is given in terms 


of the Schwarzschild radius as, 
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Here vı and r are the horizontal and vertical viscosities 
discussed earlier in §VI.A.2. 


If the conditions are right to trigger it, calculations 
show that the non-linear evolution of the instability un- 
der X-ray binary conditions leads to evolution that is 
qualitatively consistent with observations (Wijers and 
Pringle, 1999). It is not so obvious, however, whether 
the conditions are right. The hydrodynamic expectation 
that v2/vı ~ 1/(2a7) > 1 (equation 317) implies that 
accretion disks are quite “stiff’, and would only be un- 
stable at unreasonably large radii for typically assumed 
values of a and 7. For example, even assuming quite large 
values of a = 0.3 and 7 = 0.3, instability only sets in for 
r > 3x 10t rg. An interesting possibility is that anal- 
ogous instabilities might be driven by the larger torques 
from pressure or magnetic forces if there is a wind from 
the disk (Lai, 2003; Schandl and Meyer, 1994). 
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D. Other conjectured instabilities 


Many other accretion disk instabilities have been sug- 
gested, with varying degrees of theoretical and / or ob- 
servational evidence. Here are a few. 

Gravo-magneto limit cycle instability. Young Stel- 
lar Objects (YSOs) show eruptive behavior, known as 
FU Orionis outbursts, that looks very much like a slowed- 
down version of dwarf nova outbursts (Audard et al., 
2014). Thermal instability models for FUOrs are pos- 
sible (Bell and Lin, 1994), but require high disk masses 
and unusually low values of a. A modified limit cycle, 
in which non-ideal MHD processes strongly damp turbu- 
lence on the cold branch, and the onset of self-gravity 
triggers the transition to the hot branch, has been sug- 
gested for these systems (Armitage et al., 2001; Gammie, 
1999b; Martin and Lubow, 2011; Zhu et al., 2009). 

Thermal wave or irradiation instability. A disk whose 
thermal structure is dominated by irradiation from a cen- 
tral source has a scale height h(r) whose profile is fixed 
by the shape of the absorbing surface of the disk. For 
YSOs and protoplanetary disks, consistent equilibrium 
disk shapes are modestly flared to large radius (Chiang 
and Goldreich, 1997; Kenyon and Hartmann, 1987). It 
is easy to imagine how such an equilibrium might be un- 
stable. If the disk surface develops a small “ripple”, the 
inward-facing part of the perturbation will intercept more 
radiation from the central source, and the outward-facing 
part less. The ripple will strengthen—perhaps to the 
point that an actual shadow is cast—and propagate ra- 
dially. This is an old idea (Cunningham, 1976; D’ Alessio 
et al., 1999; Dullemond, 2000; Watanabe and Lin, 2008) 
that has received renewed study (Ueda et al., 2021; Wu 
and Lithwick, 2021) as a potential cause of the annular 
structures observed in protoplanetary disks. 

Magnetic Prandtl number driven instability. The mi- 
crophysical viscosity v and resistivity 7 of a plasma have 
the same dimensions, and from them one can construct 
a dimensionless ratio, the magnetic Prandtl number, 
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Pm (376) 


Numerical simulations, both of isotropic stirred MHD 
turbulence (Schekochihin et al., 2004) and of the MRI 
in shearing boxes (Fromang et al., 2007), exhibit a de- 
pendence on Pm. The modest value of fluid Reynolds 
number realized numerically means that these results 
must be approached cautiously, but they are consistent 
with a plausible physical argument in which the value 
of the Prandtl number affects the rate of reconnection 
and hence large-scale fluid properties. Balbus and Henri 
(2008) noted that the value of Pm in X-ray binary disk 
models goes from Pm > 1 close to the black hole to 
Pm < 1 at 10? rg, which might feed forward into observ- 
able time-dependent behavior (Potter and Balbus, 2014). 
A possible confounding factor is the importance of radi- 
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ation viscosity in the relevant disk regions. 


VIII. GEOMETRY OF ACCRETION 


In several circumstances, the existence of a particular 
geometry (spherical symmetry, axisymmetry) allows us 
to calculate how gas flows from large radii toward the 
accreting object. We reproduce here a couple of well- 
known results, focusing on the simplest cases. 


A. Spherical (Bondi) accretion 


A point mass M is at rest with respect to gas, uni- 
form at infinity, with density pə and sound speed cs. 
The problem is spherically symmetric and the mass ac- 
cretion rate can be estimated from dimensional analysis. 
A spatial scale can be constructed as, 


GM 


= 
c3 


Rg = (377) 
called the Bondi radius. Multiplying the spherical sur- 
face area 47 R2 by the characteristic density and by the 
characteristic velocity, the accretion rate is, 


GM) 
( 1) Ta 


Mp ~ 4r (378) 
This is the Bondi accretion rate. 

The dimensional analysis gives the correct estimate, 
but we can solve the problem exactly (Bondi, 1952). We 
need the continuity and momentum equations, which for 
a steady spherically symmetric flow read, 


M = Arr Upp, 


dv idp GM 
= 
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We also need an equation of state. The simplest case 


is the isothermal one, with p = pc?. The momentum 
equation is then, 


2d mer 2 dnp GM 

T dr ° dr r2 ` 

We take the log of the continuity equation and differen- 
tiate with respect to radius to obtain, 


(381) 


2 dln dlnv 
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r dr dr 
and use this relation to eliminate density from the mo- 
mentum equation. The result is, 


(v2 — e) dlnv, 2c (1 =) . 


Cc = 
nS! dr r 2rc2 
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This equation can describe a number of physical circum- 
stances, but the one that is almost always of interest as- 
trophysically is the one where the flow speed —> 0 at large 


radius, and approaches some large value >> c, near the 
central object (this must be true for a black hole, and is 
usually true even for objects with surfaces). There must 
therefore be a radius, called the sonic radius rs, where 
v2 = c2, where the flow makes a transition between sub- 
sonic and supersonic motion. Inspection of the above 


equation pins the sonic radius at, 


GM 


Ts = =. 
2 
2c? 


(384) 
Writing ps = p(rs), the accretion rate for an isothermal 
equation of state is just M = 4rr2c.ps. 

To complete the solution, we need to relate the density 
at the sonic point ps to the density at infinity p.. This 
is accomplished using the Bernoulli constant, 

1 dp GM 


H = 2v? . 
ur + 7 7 


(385) 


For an isothermal gas at temperature T we have that 
dp = (R/u)Tdp, so H = (1/2)v2 + cln p — GM/r. The 
relation between the general conditions in the flow and 
those at the sonic point is then, 


1 GM GM 
v+ elp =- + emps- . (386) 
2 r 2 Ts 
On substituting for rs, 
2GM 
v = 22 [m (2) J pei, (387) 
p 2 r 


Noting that as r — oo the radial velocity vanishes, we 
find that In(ps/po.) = 3/2 and the accretion rate is, 


(GM)? 
3 


M = re’? Po: (388) 


S 


This differs by about 10% from the answer suggested by 
dimensional analysis. 

The extension to an adiabatic equation of state is not 
hard, and is covered in many textbooks (e.g. Clarke and 
Carswell, 2007; Frank et al., 2002). The same analysis 
applies equally to the hydrodynamics of inflows and out- 
flows, and in the latter case is the basis of the Parker 
wind solution that describes a thermally driven stellar 
wind (Parker, 1958). 


B. Supersonic accretion from an external medium 


A second limit amenable to analytic analysis occurs 
when a point-mass accretor moves supersonically through 
a medium that is uniform at infinity. Examples include 
an isolated neutron star or black hole traversing a molec- 
ular cloud, or a compact object accreting from the wind 
of a massive star in a binary system. This is Bondi- 
Hoyle-Lyttleton accretion (Bondi and Hoyle, 1944; Hoyle 
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accreted <e —> not accreted 


accretion 
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cancellation 
of non-radial 
momentum 


ballistic trajectory 


FIG. 25 Illustration of the physics of Bondi-Hoyle-Lyttleton 
accretion. A point mass moving supersonically through a uni- 
form medium focuses streamlines of the fluid toward a trailing 
accretion line. Streamlines collide in diametrically opposed 
pairs when they reach the accretion line, retaining only their 
radial momentum. Close to the point mass, the total energy 
after collision is negative, and the gas accretes. 


and Lyttleton, 1939). Although it is, in principle, at least 
an axisymmetric fluid problem, a calculation that largely 
ignores pressure effects is instructive and physically use- 
ful. Edgar (2004) provides a clear review, and we follow 
his treatment of the basics here. 


The physical argument for how Bondi-Hoyle-Lyttleton 
accretion proceeds is illustrated in Figure 25. The flow 
is evidently axisymmetric at large distances from the ac- 
creting object, and we assume that there are no insta- 
bilities that spoil this symmetry closer in. We further 
assume that, because the motion is supersonic, fluid fol- 
lows essentially ballistic trajectories in the gravitational 
field of the point mass. Given these assumptions, pairs of 
streamlines that pass by the mass on opposite sides are 
focused toward collision points that lie on a straight line 
behind the accreting object. The collision cancels out the 
non-radial component of the momentum, leaving only a 
radial component v(r) (which we will calculate shortly). 
The specific energy, (1/2)v? — GM/r, is an increasing 
function of distance r along the accretion line. Gas with 
negative specific energy after collision flows back along 
the accretion line and gets accreted, while that further 
away flows outward. The accretion rate is then given by 
the mass flux into a cylinder, whose size is defined by the 
impact parameter of the critical streamline that has zero 
specific energy after collision on the accretion line. 


Quantifying the above argument requires calculating 
the impact parameter of the hyperbolic orbit that leads 
to zero energy after collision along the accretion line 
(Edgar, 2004). The setup for the calculation is shown 
in Figure 26. A point mass M moves with velocity vy 
with respect to a uniform medium of density px. We 
work in polar co-ordinates (r,@), and consider the orbit 


impact 
parameter b 


FIG. 26 Setup for the calculation of the Bondi-Hoyle- 
Lyttleton accretion rate. 


of a point particle with impact parameter b. Ignoring 
fluid effects, the equations of motion are, 


ee GM 
T= ro = ~ ge a (389) 
7? = bvo. (390) 


The dots denote time derivatives. The specific angular 
momentum h = bua. 

Solving the above equations to determine the shape of 
the orbit is a standard exercise in celestial mechanics. 
Substituting u = r—', and repeatedly using the chain 
rule, one obtains, 


—~+u=— >. 391 
dé? h? a) 
The solution of this ODE can be written in various ways. 
Often one works toward the real space solution, 
h?/GM 
r= a (392) 
1+ ecos (0 — w) 
describing a general conic section, with e and w being 
constants that one identifies with the orbital eccentricity 
and longitude of pericenter. For our purposes, it’s a bit 
simpler to instead write the solution as, 


. GM 
u = cı cos 0 + c2 sin 0 + ——, 


3 (393) 


with cı and cz being constants. The constants are deter- 
mined by requiring that as 0 > 7, 


u —> 0, (394) 

i vo. (395) 
Noting that, 

paap" (396) 
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the solution is, 


M 
u= [1 + cos 6] — > sing, 
Finally, setting 0 = 0 we find that the streamlines with 
impact parameter b collide along the accretion line at a 
distance, 


(397) 


Pan bv 

2GM 

downstream behind the point mass. Post-collision, the 

azimuthal velocity goes to zero, leaving only a radial com- 
ponent, 


(398) 


(399) 


P= Uoo: 


Note that the radial velocity is independent of the dis- 
tance of the collision point from the point mass. 

Using the above results, we assess whether the post- 
collision gas is bound to the accretor. The condition for 
being bound is that, 


(400) 


The critical streamline that yields just-bound gas has an 
impact parameter, 
2GM 


berit = 2° 
Uso 


(401) 


All gas flowing into the accretion cylinder with b < berit 
will be accreted, so the accretion rate is, 
(GM)? 


3 Pœ: 
co 


M = 12 4 VooPoo = 4m (402) 
Up to a numerical factor that depends on the adiabatic 
index, the Bondi-Hoyle-Lyttelton accretion rate is just 
the Bondi accretion rate, with the relative velocity of the 
moving accretor replacing the sound speed that enters 
into the spherical formula. 

Numerical simulations of Bondi-Hoyle-Lyttleton accre- 
tion, which is a computationally hard problem, go back 
a long way (Livio et al., 1986; Matsuda et al., 1991; Ruf- 
fert, 1999; Taam and Fryxell, 1988). The analytic es- 
timate for the accretion rate is found to be a reason- 
able estimate, but the flow morphology depends on the 
adiabatic index, features a bow shock, and is often un- 
steady. Going beyond the simplest version of the prob- 
lem, the rates of mass and angular momentum accre- 
tion from media with density or velocity gradients, or 
stochastic inhomogeneities, are of interest. These more 
complex situations are hard to treat analytically (for a 
discussion of the physical considerations, see Davies and 
Pringle, 1980). Relatively recent works include those 
by Krumholz et al. (2006), Blondin and Raymer (2012), 
MacLeod and Ramirez-Ruiz (2015), and Xu and Stone 
(2019). 
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